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Abstract 


The behavior of the velocity gradient tensor, A i} = dui/dij, was studied using three 
turbulent flows obtained from Direct Numerical Simulation. The flows studied were: 
an inviscid calculation of the interaction between two vortex tubes, a homogeneous 
isotropic flow, and a temporally evolving planar wake. Self-similar behavior for each 
flow was obtained when was normalized with the mean strain rate. The case of 
the interaction between two vortex tubes revealed a finite-sized coherent structure 
with topological characteristics predictable by a Restricted Euler model. This struc- 
ture was found to evolve with the peak vorticity as the flow approached singularity. 
Invariants of A^ within this structure followed a straight line relationship of the form: 
A 3 + + r _ q where Q and R are the second and third invariants of A l3 , and the 

eigenvalue A is nearly constant over the volume of this structure. Data within this 
structure have local strain topology of unstable-node/saddle/saddle. The character- 
istics of the velocity gradient tensor and the anisotropic part of a related acceleration 
gradient tensor H tJ were also studied for a homogeneous isotropic flow and a tempo- 
rally evolving planar wake. It was found that the intermediate principal eigenvalue 
of the rate-of-strain tensor of Hij tended to be negative, with local strain topology of 
the type stable-node/saddle/saddle. There was also a preferential alignment between 
the equivalent vorticity vector and this intermediate principal eigenvalue direction. 
The magnitude of H tj in the wake flow was found to be very small when data were 
conditioned at high local dissipation regions. This result was not observed in the rel- 
atively low Reynolds number simulation of homogeneous isotropic flow. A Restricted 
Euler model of the evolution of A i} was found to reproduce many of the topological 
features identified in the simulations. 
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Chapter 1 


Introduction 


1.1 Background 

Turbulence is one of the most complicated types of fluid motion. Turbulent flow is 
generally random, disordered and unpredictable, yet it is governed by the Navier- 
Stokes equations. Almost every fluid motion in nature, engineering applications, 
and everyday life contains turbulence. The magnificent photosphere of the Sun, the 
gigantic cloud systems that affect the Earth’s weather pattern, the boundary layers 
growing on aircraft wings, the wakes of submarines, the pouring of cream into a cup 
of coffee, and the flight of a golf ball; all may contain turbulent motions. 

In 1937, G. I. Taylor and T. Von Karman[16] defined turbulence in the following 
way: 


Turbulence is an irregular motion which in general makes its appearance 
in fluids, gaseous or liquid, when they flow past solid surfaces or even 
when neighboring streams of the same fluid flow past or over one another. 

Hinze[14] formulated a more precise definition of turbulence: 

Turbulent fluid motion is an irregular condition of flow in which the various 
quantities show a random variation with time and space coordinates, so 
that statistically distinct average values can be discerned. 
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According to Hinze, it is not sufficient to define turbulent motion as being irregular 
in time alone, but in space as well. Therefore, it is impossible to describe turbulent 
motion in all detail as a function of time and space coordinates. Instead, one has to 
rely on laws of probability and statistics. It is important to note that for any study of 
turbulence to be possible, one has to accept that no matter how confusing turbulence 
is, it must still be governed by laws of physics. 


1.2 Time- averaged turbulence models 

In the time-averaged (Revnolds-averaged) description of turbulence, the number of 
unknown variables is larger than the number of equations. At present, most turbulent 
models rely on intuition, experience and dimensional analysis in an attempt to close 
the system of governing equations. 

Kline, Cantwell & Lilley[19] concluded that despite intensive computer-assisted 
activity in turbulence research, none of the existing modeling methods could predict 
a worthwhile range of perturbed shear layers with a constant set of empirical co- 
efficients. Finding a “universal” turbulence model capable of predicting any flow to 
acceptable engineering accuracy remains a challenge to research scientists world- wide. 
Ferziger, Kline, Avva, Bordalo & Tzuoo[10] modified Kline’s “zonal” modeling into a 
framework for adjusting the coefficients of a chosen turbulence model from one part 
of the flow to another. A zone is defined as a region in which a chosen turbulence 
model will perform acceptably with a given set of coefficients, which may be func- 
tions of local parameters rather than absolute constants. It is usually an identifiable 
species of turbulent flow (e.g. boundary layer, wake or jet) or a strong perturbation 
of a given species, such as a shock- wave/boundary interaction. Zonal modeling is 
needed because current turbulent models cannot accurately simulate all the flows 
with a single set of coefficients. One of the most widely used turbulent models today 
is the “k - e” model[26, 27]. Using transport equations for both the turbulent kinetic 
energy k and dissipation s, this model predicts some aspects of the behavior of tur- 
bulent flows adequately within its assumptions. Many existing turbulent models are 
variations of the “k - e” model. The purpose of these models is to solve the Reynolds 
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averaged Navier-Stokes equations in the hope of predicting the general behavior of 
turbulent flows. This behavior depends strongly on the turbulent kinetic energy and 
dissipation, the later being governed by the instantaneous velocity gradient dui/dxj 
in the flow. A good understanding of this instantaneous velocity gradient tensor is 
therefore essential to the studies of turbulence. 


1.3 Motivation 


Instantaneous velocity gradient tensor, defined as A+j = dui/dxj , in a turbulent flow 
governs the behavior of the turbulent kinetic energy dissipation, e. The transport 
equation for the turbulent kinetic energy is: 


Dq 2 

~Dt 


_d_ 

dx k 


u k ( P + Q 2 ) ~UiU k 


dU, 

dx k 




d 2 u 2 i d 2 UiU k 

dxidx k 




+ 




diij 

dxi 


(i.i) 

q 2 = UiUi /2 is the turbulent kinetic energy while V, V, V and e are the convective 
diffusion, production, viscous diffusion and mean rate of dissipation of turbulent 
kinetic energy respectively. 

Among these terms, e is of the greatest interest since it is governed by the velocity 
gradient tensor, A tJ -. 


dUi , dUi &Uj . j—. — r- 

£ = u d^ { d^ + d^ ) = " Aij{Aij + Aji) - 


( 1 . 2 ) 


Dimensional analysis is applied to the transport equation for turbulent kinetic 
energy to obtain the relative scaling of Aij. Both velocity fluctuation u,, and the 
mean velocity Ui, scale with the free stream velocity U 0 . Using this velocity scale U Q 
together with the mean flow length scale S, the relative scaling of production V is: 

dUi Ul 


V = UiUfc- . 

ox k o 

In most homogeneous shear flows, V ~ £. Hence, 

U z 

E — V Aij^Aij “f Aji) ^ ^ * 


(1,3) 


(1.4) 
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This relation implies that: 



UoS ui 


£2 R e& 



(1.5) 


Therefore, 

(1.6) 

Dimensional arguments applied to the transport equation for turbulent kinetic energy 
indicates that the instantaneous velocity gradient is larger than the mean gradient 
U a /5 by at least a factor of Re 1 / 2 . 

In addition to dissipation, velocity gradient tensor also governs the behavior of 
vorticity oj, in turbulent flows: 


^ ijkAf-j , 


(1.7) 


where e ijk is the alternating unit tensor such that: 


{ +1; if ijk are in cyclic order 123123 

— 1; if ijk are in anticyclic order 321321 (1.8) 

0; otherwise 

Figure 1.1 shows a contour plot of iso-enstrophy (magnitude of vorticity) surfaces 
in a direct numerical simulation of an evolving planar wake[30). Free stream flows 
from left to right, as indicated by the arrow. A cross section (y-z plane) is cut 
across the x-axis to illustrate the microscale regions. Figure 1.2 shows the contour 
plots of enstrophy and dissipation along this cross section, looking along the flow 
direction. The highest contour is indicated in bright red. It is observed that the 
highest intensities of both vorticity and dissipation (hence velocity gradients) occur 
in regions with very small length scale compared to the mean flow. 


1.4 Previous studies of the velocity gradient ten- 
sor 

A considerable amount of research has been directed at understanding the instan- 
taneous velocity gradient in turbulent flows. Large instantaneous velocity gradients 
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Figure 1.1: Contour plot of iso-enstrophy surfaces in an evolving plane wake. — ► 
indicates free stream flow direction. 
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(b) 


Figure 1.2: y-z plane in the evolving planar wake showing contour plots of (a) local 
enstrophy density, (b) dissipation. 
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usually occur in small-scale structures whose characteristic length scales are much 
smaller than the mean flow, as shown in figure 1.2. 

Recent studies of small-scale structures using direct numerical simulations (DNS) 
and laboratory experiments suggested common behavior in the small scale velocity 
gradients; even when the large-scale motions are very different. Ashurst, Kerstein, 
Kerr & Gibson[l] studied the behavior of the velocity gradient in an incompressible, 
isotropic homogeneous flow in a periodic box by analyzing the velocity gradient ten- 
sor (each tensor consists of nine components of dui/dxj ) at each grid point in the 
flow. The study was done using direct numerical simulation, with Taylor microscale 
Reynolds number Re x ~ 83. By looking at the probability distribution function of the 
cosine of the angle between the vorticity vector and the three principal eigenvectors 
of the rate-of-strain tensor (the symmetric part of the velocity gradient tensor), they 
concluded that there was a strong preferential alignment of the vorticity vector with 
the intermediate eigenvector. This same conclusion was also reached by studying an 
incompressible homogeneous shear flow with similar Reynolds number. They also 
found that the principal eigenvalues of the rate-of-strain tensor (sorted in descending 
order such that a > (3 > 7) have a preferred ratio ofa:^: 7«3 :l: — 4 in regions 
of high dissipations. 

Vincent & Meneguzzi[35] performed a direct numerical simulation of a homo- 
geneous isotropic turbulence at a higher Reynolds number, with Re x ~ 150 and 
ReL ~ 1000. Visualization of the flow field at different times showed that the strongest 
vorticity in the flow was organized in very elongated thin tubes, in agreement with 
the findings of Ruetsch & Maxey[29] in a similar flow with Re x « 60. Vincent k 
Meneguzzi also concluded from this study that there is a preferential alignment of 
the vorticity vector with the intermediate principal eigenvector of the rate-of-strain 
tensor. 

Soria, Sondergaard, Cantwell, Chong k Perry[31] studied the dissipating motions 
of an incompressible mixing layer with both laminar and turbulent initial conditions, 
at Res ~ 250. Their results showed that regardless of initial conditions, the bulk of 
the total kinetic energy dissipation is contributed by fluid structures with local strain 
rate topology characterized as unstable-node/saddle/saddle (a > {3 > 7; /? > 0). 



8 


CHAPTER 1. INTRODUCTION 


Blackburn, Mansour & Cantwell [2] investigated the topological features of the 
velocity gradient field in turbulent channel flow. The Reynolds number of this simu- 
lation was Re = 7860, based on half channel width and the center-line velocity. In all 
regions of the flow, there was a strong preference for the vorticity to align with the 
intermediate principal eigenvector of the rate-of-strain tensor. Away from the wall 
regions, the intermediate principal eigenvalue of the rate-of-strain tensor (these prin- 
cipal eigenvalues are also known as strain rates) tends to be positive at sites of high 
viscous dissipation of kinetic energy. The velocity gradient tensor also showed pref- 
erence for local flow topologies of unstable-node/saddle/saddle. They also used the 
discriminant of the velocity gradient tensor to identify flow structures which extends 
from very close to the wall to the free stream. 

Tsinober, Kitt & Dracos[32] conducted laboratory experiments on turbulent grid 
flows and on the turbulent boundary layer over a smooth plate. To determine the 
invariant properties of the flows, all nine components of the velocity gradient tensor 
at every grid point in the flow were measured using a hot wire probe. In both flows, 
the probability density distribution of the alignment angle between the vorticity and 
the three principal eigenvectors of the rate-of-strain tensor was determined. As in 
the DNS studies mentioned above, they confirmed that there was a strong tendency 
for the vorticity to align itself with the intermediate eigenvector of the rate-of-strain 
tensor. 

The DNS and experimental studies mentioned above demonstrated that small- 
scale motions exhibit common behavior even when the large-scale motions of various 
homogeneous turbulent flows are different. An analytical model capable of predicting 
the behavior of these small-scale motions would be very useful. To find such a model, 
an understanding of the behavior of the velocity gradient tensor in turbulent flows is 
essential. Jimenez[15] suggested a model which described the vorticity distribution of 
a stretched vortex tube. Using this model, Jimenez was able to explain the observed 
preferential alignment of the vorticity vector with the positive intermediate principal 
eigenvector of the strain-rate tensor using purely kinematic arguments. Girimaji & 
Pope[ll] modelled the velocity gradient as a diffusion process. Using a stochastic 
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model by following fluid particles in incompressible, homogeneous isotropic turbu- 
lence, they reproduced one-time and two-time statistics in good agreement with those 
obtained from full turbulence simulations. Cantwell[4] derived an analytical solution 
to the Restricted Euler equation, first studied by Vieillefosse[33, 34]. The asymptotic 
solution reproduced the observations found in direct numerical simulations of turbu- 
lent flows: two of the principal eigenvalues are positive, and the vorticity vector aligns 
with the intermediate principal eigenvector exactly. Girimaji & Speziale[12] analyzed 
a modified version of the Restricted Euler equation numerically by including the ef- 
fects of the mean velocity gradient tensor. This model preserved the balance of mean 
momentum for most homogeneous turbulent flows with mean velocity gradients. 


1.5 Classification of local flow topology 

Many studies of turbulent flows have been performed recently. Vast quantities of data 
have been collected either using sophisticated laboratory instruments or generated by 
computer simulations. With so much information available, the need for a systematic 
and efficient way to analyze these data is crucial. Classification of local flow topol- 
ogy in turbulent flow, using critical point theory as described by Chong, Perry & 
Cantwell[8] is one of the best ways to analyze the large amounts of data associated 
with these studies. A critical point is a point in the flow field where all three velocity 
components are zero and the streamline slope is indeterminate. The flow topology at 
each point in the flow may be accessed from the viewpoint of an observer travelling 
with the local velocity of the flow. For such an observer, each point is a critical point, 
and the topology can be categorized accordingly. In addition to Chong et al., Soria et 
a/. [31], Sondergaard[30], Chen et aZ.[7], Blackburnet aZ.[2] and others have also given 
insightful descriptions of the various classifications of local flow topologies in different 
types of turbulent flows. 

Most of the results in this thesis will be expressed in terms of local flow topolo- 
gies. Therefore, the basic definitions and physical interpretations of the various flow 
topologies will be briefly described in the following sections. 
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1.5.1 Invariants of the velocity gradient tensor 

A vector field u = ui(xi,X 2 ,xz),U 2 {xi,X 2 ,x 3 ),U 3 (xi,X 2 ,x 3 ), is linearized about a 
critical point O, i.e., u(O) = 0. (x' l ,x' 2 ,x' 3 ) are the coordinates measured relative to 
the point O. Assuming that a Taylor series expansion can be carried out about the 
point O 

On ^12 °13 x \ 

u — A • x> + ■ • • = 021 0,22 O 23 x 2 + higher order terms. (1.9) 

031 032 O33 _ x ' 3 

d= A v= KT ■ (110) 

UX J Q 

Since u is the velocity field, A tJ is the velocity gradient tensor at the point O. The 
three eigenvalues, A, of A are obtained as solutions of the characteristic equation 

\ 3 + P\ 2 + Q\ + R = 0. (1.11) 

P, Q and R are the tensor invariants of A. 

1.5.2 General three-dimensional flows 

The definitions of these tensor invariants for general three-dimensional flows are given 

by: 


P = —An = —trace[A]. 

(1.12) 

Q = ^(P 2 - trace[A 2 }). 

(1.13) 

R = i(— P 3 + 3 PQ — trace[A 3 ]) = — det[A]. 

O 

(1.14) 


The characteristic equation has three different combinations of roots, (1) all real 
roots which are distinct, (2) all real roots where at least two of them are equal, or (3) 
one real root and a conjugate pair of complex roots. 

In P-Q-R space, the surface which divides the real solutions from the complex 
solutions is given by 


27 R 2 + (4 P 3 - 18 PQ)R + (4 Q 3 - P 2 Q 2 ) = 0, 


(1.15) 
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where terms have been gathered to form a quadratic equation in R. 

Chong et al. [8] gave very detailed descriptions of different types of critical points 
depending on the location of the invariants in the P-Q-R space. 

1.5.3 Incompressible three-dimensional flows 

The definitions of the three invariants for incompressible flows are much simpler due 
to the incompressibility condition, which forces the trace of the velocity gradient 
tensor, A, to zero. This condition greatly simplifies the expressions of the various 
invariants: 

P = —An = —trace[A\ =0. (1-16) 

Q = “(trace [A 2 ]) = -^A im A mi . (1.17) 

R = ~^{trace[A 3 ]) = -det[A] = "AimAmkAu. (1.18) 

Because the first invariant P is now zero for incompressible flows, the local flow 
geometry is completely determined by the location of the second and third invariants 
of the velocity gradient in the two-dimensional Q-R space. 

The surface that separates real eigenvalues from complex eigenvalues now reduces 
to a curve given by: 

D = (27/4)i? 2 + Q 3 = 0. (1.19) 

The value of the discriminant, D determines the nature of the eigenvalues of A. When 
D > 0, A has one real and two complex eigenvalues. D < 0 implies that A has three 
distinct, real eigenvalues. When R = ±(2-\/3/9)(— Q) 3 ^ 2 and D = 0, A has three 
real eigenvalues, two of which are equal. The curve D = 0 acts as a dividing line 
which separates solutions with three real eigenvalues from those with one real and 
two complex conjugate eigenvalues. 

The curve D = 0 and the Q axis separate the Q-R space into four regions. Fig- 
ure 1.3 illustrates the different topologies in these four regions: 

1. above the dividing curve and to the left of Q axis, the local flow spirals in 
towards the local origin in a plane and then flows out along the third direction. 
This local flow geometry is referred to as stable-focus/stretching. 
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Figure 1.3: Three dimensional topologies in Q-R ( P = 0) invariant space 

2. above the dividing curve and to the right of Q axis, the local flow approaches 
the origin along one axis and spirals out in a plane. This local flow geometry is 
referred to as unstable-focus/compressing. 

3. below the dividing curve and to the left of Q axis, the local flow approaches 
the origin along two axes and flows outward along the third. This local flow 
geometry is referred to as stable-node/saddle/saddle. 

4. below the separator and to the right of Q axis, the local flow approaches the 
origin along one axis and flows outwards along the other two. This local flow 
geometry is referred to as unstable-node/saddle/saddle. 

1.5.4 Construction of a typical Q-R invariant plot 

Figure 1.4 illustrates a typical Q-R plot obtained for a direct numerical simulation 
of incompressible, homogeneous isotropic flow. At one instant in the simulation, the 
values of Q and R at each grid point are calculated using the definitions given in the 
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Figure 1.4: A typical contour plot of joint pdf of Q vs R. 


previous section and cross-plotted to the Q-R space. This process is repeated for all 
grid points in the whole flow. Figure 1.4 plots the contours of the number density of 
points lying within a unit area in the invariant space. Logarithmic contour levels (1, 
10, 100 etc.) are chosen due to extreme variations of the number density that occurs. 
In this way, isolated points far from the origin corresponding to the largest gradients 
in the flow (which occur in the smallest-scale structures) are captured with contour 
level of 1. At the same time, information about the distribution of invariants near the 
origin corresponding to intermediate-scale and large-scale structures is also provided 
by the higher contour levels. This invariant plot becomes the joint probability density 
function (pdf) between Q and R when the contour levels are normalized by the total 
number of data points in the plot. 
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1.5.5 Invariants of the rate-of-strain and rate-of-rotation ten- 
sors 

Any velocity gradient tensor A t j may be decomposed into the rate-of strain tensor S 
and the rate-of-rotation tensor W. The components of S and W are defined as : 


$ij — {Aij + Aji)j 2. 

(1.20) 

W{j = ( Aij — Aji)/2 . 

(1.21) 


Considering only incompressible flows, the invariants for S and W are : 


P s = -Sa = -trace[S] = 0. (1.22) 

Qs = "•S'irnS'rm = ~trace[S 2 ]. (1.23) 

Rs = “ SimSmkSki - ~trace[S 3 ]. (1.24) 

P w = —Wa = — trace[W ] = 0. (1.25) 

Qw = -^W im W mi = -^trace[W 2 ]. (1.26) 

Rw = -^W im W mk W ki = 0. (1.27) 

The invariants of the velocity gradient tensor are related to the invariants of the 
rate-of-strain and rate-of-rotation tensors in the following ways: 

Q = Qs + Qw (1.28) 

R = Rs- W im w mk s ki . (1.29) 


Since 5 is a symmetric tensor, all the eigenvalues of S must be real. Therefore, 
all the data points must lie below the discriminant curve D s = (27/4 )R 2 + Q* = 0 in 
the Q s -R s invariant space. 
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1.5.6 Physical meanings of various invariants 

The invariants obtained from the A, S and W tensors have direct implications for the 
fluid flows. In addition to classifying the local flow fields as observed by an observer 
travelling with the fluid, these invariant plots also reveal physical information about 
the flow: 

1. the second invariant of S, Q s , is directly proportional to the local dissipation. 

e = 2 i/SijSji = —4 vQ s . (1.30) 

Points with large negative values of Q s correspond to high dissipation regions 
in the flow. 

2. the second invariant of W, Q w , is the magnitude of the vorticity squared, also 
known as the local enstrophy density. 

Qw = mW mi = ViUJi. (1.31) 

Points with large Q w have high enstrophy density. 

3. the difference between the third invariant of the rate-of-strain tensor and the 
third invariant of the velocity gradient tensor, ( R s — R), reveals information 
about the vortex stretching rate, a. 

(R s ~R) = W im W mk S ki = a. (1.32) 

Information regarding the physical flows can therefore be obtained by inspecting 
the invariant plots. If Q is large and positive, the local enstrophy density is large 
and dominates the local dissipation ( Q w » Q s ), as in the case of solid body rotation 
near the center of a vortex tube. On the other hand, when Q is large and negative, 
the local dissipation is large and dominates the local enstrophy density ( Q s 2> Q w ), 
as found in flows with pure straining. Finally, when (Q s ~ Q w ), the local dissipation 
and local enstrophy density are comparable in magnitude, as found in vortex sheets. 
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1.6 Outline of thesis 

The main objective of this thesis is to study the characteristics of the velocity and 
acceleration gradient tensors in one inviscid flow simulation and two viscous flow 
simulations. The inviscid flow is chosen because it exhibits behavior consistent with 
the Restricted Euler model of the evolution of Ay. All the flows studied in this thesis 
are obtained from direct numerical simulations, where the Navier-Stokes equation is 
solved exactly using “spectral” methods. There is no attempt to model any unresolved 
scales in these simulations. 

Chapter 2 describes the Restricted Euler equation and its analytical solution. 
The definition of the acceleration gradient tensor, which governs the behavior of 
the velocity gradient tensor, is presented. A velocity gradient tensor with random 
components obtained from a Gaussian distribution with zero mean and unit variance 
evolving with the Restricted Euler model is also discussed. 

Chapter 3 analyses the local flow topologies of an incompressible Euler calculation 
of two interacting vortex tubes simulated by Kerr[17]. Figure 1.5(a) depicts the 
domain surrounding the peak vorticity in this flow. Different regions with interesting 
local flow topologies are identified. 

Chapter 4 describes the numerical approach used by the author to perform a direct 
numerical simulation of a homogeneous isotropic flow. A contour plot of iso-vorticity 
surfaces of this simulation is shown in figure 1.5(b). The behavior of the velocity and 
acceleration gradient tensors in this flow is analyzed using classification of local flow 
topologies. 

Chapter 5 presents the results obtained from a simulation of a temporally evolving 
plane wake generated by Sondergaard[30]. Behavior of the velocity and acceleration 
gradient tensors in this relatively higher Reynolds number simulation is compared 
and contrast with results obtained from the homogeneous isotropic flow. A general 
view of this flow has been shown in figure 1.1. 

Chapter 6 gives the major conclusions and outlines some recommendations for 
future studies. 
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Chapter 2 

Equations for the velocity gradient 
tensor 


2.1 Introduction 


Ashurst et a/.[l] examined the statistical properties of the velocity gradient tensor 
in a forced isotropic turbulence and in a homogeneous shear flow, using DNS. They 
observed that there is a strong tendency for the velocity gradient tensor to approach 
a state where two of its principal rates-of-strain are positive while the remaining one 
is negative. There was also a preferential alignment between the vorticity vector and 
the intermediate principal eigenvector of the rate-of-strain tensor. Both tendencies 
became more pronounced when data were conditioned at higher dissipation regions. 
Similar observations were also reported in other studies of homogeneous turbulent 
flows. 

To understand the behavior of velocity gradient tensor Aij = dui/dxj, the trans- 
port equation for A y needs to be studied closely. The derivation of this transport 
equation is described below. 

The incompressible Navier-Stokes equations: 
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are differentiated with respect to Xj leading to 


dAjj Aij 

H Uk ~ — h A ik A k j = 


d 2 p d 2 A^ 

+ v- 


( 2 . 2 ) 


dt dx k ^ dxidxj ' dxkdxk 

For incompressible flow, An = 0. Therefore, the pressure is given by taking the 


trace of the above equation: 


AikAki — 


d 2 p 


dx k dx k 

Subtracting equation 2.3 from equation 2.2 yields 


dAij dA.a 
+ u k 


” + AikAkj — AkmAmk-^- = H, 


dt dxk 

where 5y is the Kronecker delta, and 

tt ( &P _ 9 2 p Sjj 

v l dxidxi dxkdxk 3 


v 


+ v 


&Aj 

dx k dxk ' 


(2.3) 


(2.4) 


(2.5) 


The Hij tensor is composed of the cross derivatives of the pressure field and the 
viscous diffusion of the velocity gradient tensor. This tensor shall be termed the 
“acceleration gradient tensor” in the rest of this thesis, even though it is not the true 
acceleration gradient tensor. A true acceleration gradient tensor, H^ can be derived 
by expanding a vector field at a critical point using a Taylor series, the same way as 
the velocity gradient tensor Ay was derived in equation 1.10. 



( 2 . 6 ) 


Differentiating this equation with time t, and using equation 2.1: 

cPxi dui dp d 2 Ui 

dt 2 dt dxi^ V dxkdxk a * 

Therefore, true acceleration gradient tensor 

deg _fr _ d 2 p d 2 A ij 

dxj ,J dx^Xj V dxkdxk 


(2.7) 


( 2 . 8 ) 


Close inspection of the true acceleration gradient tensor reveals that Hy is almost 
identical to Hy , except that Hy has been forced to be trace-free, which makes the 
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analysis of this tensor easier. Therefore, H v is essentially the an-isotropic part of the 
true acceleration gradient tensor. The behavior of this Hij tensor will be discussed 
in later chapters. 

The evolution equations for the invariants of the velocity gradient tensor are de- 
rived by forming the appropriate double and triple products with equation 2.4 and 
taking the trace: 

^ + 3 R=-A ik H ki , (2.9) 

at 

d 4-\Q 2 = -AinA nm H mi . ( 2 . 10 ) 

at o 

2.2 The Restricted Euler model 

Cantwell [4] analyzed the homogeneous case with Hij = 0, which will be referred 
to as the Restricted Euler model for Aij. This equation was first studied by Vieille- 
fosse [33, 34]. Setting H t j — 0 converts equation 2.4 into a system of coupled ordinary 
differential equations (ODE) for all nine components of Aij(t). Under this assump- 
tion, the velocity gradient tensor of a fluid particle evolves independently of other 
particles in the flow since the effects of the pressure and viscosity have been removed 
from the governing equation. This system of ODE’s was solved analytically. A brief 
description of the derivation of this analytical solution is presented below. 

Given an initial velocity gradient tensor, A ij, its discriminant D is related to its 
second and third invariants such that: 

D = Q 3 + ^R 2 = Q\ + j Rl = Ql (2.11) 

The quantities Qi and are the initial values of Q and R respectively. Q 0 is the 
value of Q when R = 0 and it is used to define the time scale for the evolution of 

Aij{t) 

t 0 =\ V abs< R°y Q° ^ 0 (2.12) 

l Qo = 0 

All relevant variables are non-dimensionalised by t a : 


Q — f — Rtcn j — Aijt 0 , T t/to 


(2.13) 
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and equation 2.11 becomes: 

27 

— r 2 + q 3 = sgn(Q 0 ) (2.14) 

where 

f +i; Qo> o 

sgn(Qo) =| -1; Q 0 < o (2.15) 

l 0; Q 0 = 0 


2.2.1 Analytical solutions of q(r) and r(r) 


By setting Hij — 0, the non-dimensionalised forms of the coupled ordinary differential 
equations 2.9 and 2.10 become: 


(2.16) 


dq 
dr 

dr , 

Tt = ^ ■ < 2 - 17 ) 

Solutions to equations 2.16 and 2.17 are expressed in Jacobian elliptic integrals of 
the first kind, F(4 > , k), and the cosine amplitude function, cn, where 


= — 3r. 
— -q 2 . 



ds 

Vl — k 2 sin 2 s’ 


cos((p) = cn(F). 


The solutions depend on Q a . 


(2.18) 


• Q 0 > 0 


q + (r) 


(1 - vg) - (1 + V§mj(2/3^)T] 

1 - cn[(2/3‘/ 4 )rJ ’ < T < Tm "' 


r+ ( r ) = 0 < t < t+„. 

C. = (3 1/4 /2)4f (|,sin(57r/12)) . 


• Q 0 < 0 


9 ( r ) 


-(1 + VS) + (1 - v / 3)cn[(2/3 1 / 4 )r] 

1 - cn[(2/3 1 / 4 )r] ’ 


0 < r < 


(2.19) 

( 2 . 20 ) 
( 2 . 21 ) 

( 2 . 22 ) 
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r ^ = V 0 <T<T max- 

T max = (3 1/4 /2)4 F (|> sin (^)) ' 

when Q 0 = 0, the solutions depend on the initial value of r as well. 
1. r(0) < 0, 


Q°(r) = - 


r °(r) = -- 


1 


k l + (l/V5)r y 


; 0 < r < oo. 


; 0 < r < oo. 


2. r(0) > 0, 


3\/3 \1 + (1/V5)r, 

^“■(i-dW) 1 o<t<v ' 3 ' 
r “ w - 57S (— dW) ; 0<T<V5 - 


(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 


2.2.2 Analytical solution of Oy(r) 


A second order differential equation for Oy can be obtained by differentiating equa- 
tion 2.4 with respect to t and expressing the resultant equation in non-dimensional 
form: 

-£%- + |g(r)0ij = 0. (2.29) 

Since the solution for q(r) is in terms of the Jacobian elliptic integrals, the determina- 
tion of the exact expression for Oy(r) is complicated. Fortunately, r(r) is a monotonic 
function of r. Expressing q as a function of r using the relationship in equation 2.14, 
and replacing the independent variable from r to r yields: 


4 

9 


enn(D \ — —r 2 ) ^22. 

sgn(Q 0 ) 4 r J dr2 


da 2 -i 2 

4r W + ?‘= = °- 


(2.30) 


The dependence of on r is now assumed implicitly by the variable r(r). 
The solution of equation 2.30 is expressed in the simple form: 


Oij(r) = Cijfi(r) + D t] f 2 (r). 


(2.31) 
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fi(r) is an even function of r and / 2 (r) is an odd function of r. C t7 (r) and £>„(r) are 
matrix constants of integration determined by initial conditions. 


The solutions of /i(r) and / 2 (r) are again dependent on the initial value of Q a : 

• Q 0 > 0 , 



/rW = \ [(1 + r)'/» + (1 - 

(2.32) 


ftir) = ^ [(1 + - (1 - — Tp-r) 1 / 3 

(2.33) 

• Q 0 < 0, 

/f ( r ) = (l + cos | ifan -1 [(3v / 3/2)r] | 

(2.34) 


/ 2 _ (r) = (2/\/3) ^1 + sin |^tan _1 [(3v / 3/2)r]| 

(2.35) 

• Q 0 = 0, 

/i°(r) = 2 1 / 3 [(3\/3/2)r]" 2 / 3 

(2.36) 


/ 2 (r) = (2 2 / s /3v^)[(3V3/2)r] 1 / 3 

(2.37) 

The solutions for Cij(r) and D i3 (r) are expressed in terms of f, a,., 
which are the initial values of r, a t3 and da l3 /dr respectively: 

and da /dr, 


Cij = dij[?(f)] 2 (^) f + (^(a ik d kj ) + q(r)6 ij S J f 2 (r) 

(2.38) 


Dij = -dij[q(r)} 2 (^)f + - (^(a ik d kj ) + q(f)6 i} ^ /x(f) 

(2.39) 


2.2.3 Numerical procedure in obtaining Aij(t) 

In summary, given an initial velocity gradient tensor, A ij} the analytical expression of 
this tensor evolving with the Restricted Euler model after time t can be determined 
using the following procedure: 

1. Obtain the characteristic time scale, t 0 from equation 2.12. 
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2. Non-dimensionalize the velocity gradient tensor and its invariants by t a , giving 
a»i00, q(r) and r(r). 

3. Determine r from q(r) from the inverse of equation 2.19, 2.22, 2.25 or 2.27, 
depending on s<?n(Q 0 )and r(r). 

4. Add the non-dimensionalised time increment At = t/t Q to the calculated r. 

5. Obtain the new values of q(r) and r(r) from equations 2.19 to 2.28. 

6. Determine fi(r) and / 2 (r) from equations 2.32 to 2.37, depending on the value 
of r(r). 

7. Get the derivatives of f\(r) and / 2 (r) by differentiating equations 2.32 to 2.37 
analytically. 

8. Obtain Cij and Dij from the initial values of f, q(f), fi(f), / 2 (f) and their 
derivatives. 

9. Calculate the final velocity gradient tensor from equation 2.31. 

10. If dimensional final velocity gradient tensor is desired, divide the calculated 
velocity gradient tensor by its own initial characteristic time scale, t 0 . 

2.2.4 Asymptotic solution of a^(r) 

The asymptotic form of a i3 as r — > oo can be expressed as : 

fan Mr) - MM) 1/3 - (2-40) 

r(r) is a function that becomes singular in finite time. 

Kij satisfies the algebraic equation 

K im K mj + ^3 Kij - 2 1 % = 0. (2.41) 

The invariants of K tJ can be derived using equation 2.41 and the continuity con- 
straint to be: 

P(K<i) = 0.0; Q{K„) = A. Jt( Kij ) = 1.0. 


(2.42) 



26 CHAPTER 2. EQUATIONS FOR THE VELOCITY GRADIENT TENSOR 


Notice that the second and third invariants of Kij tensor lie on the boundary of 
the dividing curve since Q(Kij) and R{K tj ) satisfy the equation: 

27 

—R 2 + Q 3 = 0. (2.43) 

The tensor can also be written as the sum of a symmetric tensor and an 
anti-symmetric tensor: 

Rij | * + | k , (2.44) 

where 

Sul* 0 0 

Sij\k = 0 S22U 0 (2-45) 

0 0 — Su|* — S 22 I* 

and 

0 — ^3 1* f2 2 |* 

Wijlk = ft 3 |* 0 -Qxl* (2.46) 

— f^l|* 0 

Substituting the components of Kjj into equation 2.41 produces nine equations, which 
are used to solve for 5y|* and f2i|*. The solutions exhibit the same topological 
characteristics seen in direct numerical simulations, namely: (1) the intermediate 
strain rate S22U must be non-negative, and (2) the vorticity vector must align with 
this intermediate strain rate. 

Since Q = Q s +Q W = -3/2 2/3 for the K l3 tensor, there exists a linear relationship 
between Q s {K i:} ) and Q w {K tJ ) such that: 

-Qs(Kij) = Q w (Kij) + ^3. (2.47) 

Other relevant relationships between various invariants of K tJ tensor are: 

Rs(Kij) = —-^-^Q s (Kij) — — . (2.48) 

and 

Qw(Kij) = -^j^(R s (Kij) — R(Kij)). (2.49) 

Figures 2.1 (a-d) depict the relationships between the different invariants of the K ti 
tensor. In figure 2.1(a), all the asymptotic solutions collapse into a single point with 
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Q(Kij) = — 3/2 2//3 and R(Kij) = 1.0 — regardless of initial condition. Figure 2.1(b) 
shows the linear relationship between Q s (Kij) and R s {Kij). Notice that since Q s (Kij) 
and R s (Kij ) also satisfy the characteristic equation: 

A 3 + XQ s (Kij) + R s (Kij) = 0; A = 2 ly/3 (2.50) 

the straight line osculates the dividing curve at a single point where Q s (Kij ) = 
— 3/2 2 / 3 and R s (Kij) = 1.0. Figures 2.1(c-d) plot the linear relationships derived 
in equations 2.47 and 2.49. 

2.3 Gradient tensor evolving with the Restricted 
Euler model 

The behavior of a velocity gradient tensor evolving with the Restricted Euler model 
was observed. The components of this velocity gradient tensor were generated with 
random numbers of Gaussian distribution with zero mean and unit variance. The 
trace of this tensor was forced to zero to satisfy the continuity constraint of an in- 
compressible flow. In a homogeneous isotropic flow, the volume integral of the second 
invariant Q is zero due to the balance of pressure gradients acting on the surface of 
a control volume: 

[ QdV= [ \v 2 PdV = f l(VP)-ndS = 0 (2.51) 

J v Jv 2 J s 2 

Therefore, the velocity gradient tensor was required to satisfy this condition as well. 

To satisfy both conditions, the velocity gradient tensor A tJ is constructed from: 
Aij = Sij + Wij, where Sij is a symmetric tensor and W t j is an anti-symmetric tensor. 
The components of both Sij and W X] tensors are obtained from the random number 
generator Xi with mean E(x) = 0 and variance a 2 = 1.0. The components of Sij and 
Wij tensors are: 

S u = Xi- + x 2 + x 3 ), 

S22 = x 2 - ^(xi +x 2 +x 3 ), 

S33 = x 3 - i(x! +x 2 +x 3 ), 
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S 12 = <$21 = X 4 , 
s 13 = S31 = Xbi 

S>2Z = S32 — 3-6 > 


W n = W 22 = W33 = 0, 

w 12 = -W 21 = yfi/3x 7 , 
w l3 = -w 3X = yi/ixg, 

W 23 = -W 32 = y4/3x 9 . (2.52) 

x\,x 2 , . . . ,xg are nine different random numbers obtained from the random number 
generator. The factor ^/4/3 is necessary to force the expected value of the second 
invariant Q of the tensor to zero. The derivation of this factor is explained in 
Appendix A. 

The analytical solution of this tensor evolving according to the Restricted Euler 
model after time t is obtained using the procedure described in the previous section. 
This process is repeated for 10 5 data points to obtain the probability density dis- 
tribution. The results are presented in the form of the local topologies, assuming a 
hypothetical “flow” where every initial component of the velocity gradient tensor has 
a Gaussian distribution slightly modified to satisfy the constraints described. 

2.3.1 Evolution of 

Figure 2.2 shows the relationship between the second and third invariants of the 
velocity gradient tensor, Q and R, after evolving for time t. The velocity gradient 
tensor and all its invariants have been non-dimensionalised by the characteristic time 
scale defined in equation 2.13. Figure 2.2(a) shows that all the data points collapse 
onto three curves, depending on the discriminant of each velocity gradient tensor. 
Positive initial discriminants were all normalized into a single value of +1, lying 
above the dividing plane which separates those tensors with all three real principal 
eigenvalues from those with complex eigenvalues. This positive discriminant of -t-1 
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curve crosses the Q-axis at a maximum of 1.0. This is the result of the normalizing 
factor, which is defined to be the value of Q at R = 0. On the other hand, all the 
negative initial discriminants were normalized into a single value of —1, thus lying 
below the dividing curve. These tensors have three real principal eigenvalues. Velocity 
gradient tensors with discriminant of exactly 0 (which are extremely rare) satisfy the 
characteristic equation exactly. Therefore, these points will lie exactly on the dividing 
curve. Data points are evenly distributed on both sides of the Q-axis since all the 
data were generated randomly. Figure 2.2(b) shows the distribution of the invariants 
after evolving for t = 2.0. There is an obvious shift of data points towards positive R 
along their respective discriminant curves, due to the monotone relationship between 
R and t. The velocity gradient tensor evolves with increasing R until the solution 
becomes singular in finite time. Figures 2.2(c) and 2.2(d) illustrate the analytical 
solutions to the Restricted Euler model after evolving for t = 5.0 and t = 10.0. Most 
of the data points have now migrated towards the lower right quadrant of the Q-R 
plot, where their solutions ultimately become infinite. 

2.3.2 Evolution towards the asymptotic solution 

The evolution of the velocity gradient tensor towards the asymptotic solution is illus- 
trated in figures 2.3 to 2.6. These figures depict the evolution of the initially random 
velocity gradient tensor shown in figure 2.2, normalized by the characteristic time 
scale and local R 1//3 . 

Figure 2.3(a) shows Q and R at initial time t = 0. All the data points collapsed 
into a single vertical line, parallel to the Q-axis and cutting the R- axis at the value 
of R = 1.0 since every tensor has been normalized by its own third invariant, R. 
However, the values of Q for this initial data set are not restricted. Hence the data 
points span along the vertical line for all values of Q. Figures 2.3(b-d) show the 
evolution of this data set after time t = 5.0, t = 10.0 and t = 20.0. In addition to 
having a unique value of R — 1.0, the data points converge slowly to the asymptotic 
value of Q = — 3/2 2 / 3 as the tensors evolve in time. This behavior was predicted by 
the asymptotic solution derived for the Kij tensor in equation 2.42. 



2.3. GRADIENT TENSOR EVOLVING WITH THE RESTRICTED EULER M0DEL31 



Figure 2.2: Time evolution of Q vs R. Aij normalized by local discriminant so that 
all data points must lie on three separate curves with discriminant D = +1, D = 0 
or D = -1. (a) t=0.0. (b) t=2.0. (c) t=5.0. (d) t=10.0. 
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The invariant plots of Q s vs R s at t = 0.0, t = 5.0, t = 10.0 and t = 20.0 are 
shown in figures 2.4. Initially, the data points distribute evenly below the dividing 
curve (the rate-of-strain tensor is symmetric and so all principal eigenvalues are real). 
However, for any velocity gradient tensor such that R s is close to R (for example, a 
symmetric velocity gradient tensor), normalizing the velocity gradient tensor by R l/Z 
forces R s to 1.0. The consequence of this appears on the plot as a subtle structure 
near the vertical line at R s = 1.0. The data points approach the asymptotic solution 
described in equation 2.48 as time increases. Figures 2.5 and Figures 2.6 show the 
time evolution of — Q s vs Q w and Q w vs (R s — R) respectively. Notice that the data 
points in figure 2.6(a) do not distribute evenly about the ( R s — R) = 0 axis. There is 
a structure along the R s - R = -1.0 vertical line because R has been normalized to 
1.0. 

2.3.3 Evolution of dimensional .Ay- 

Figures 2.7 to 2.10 show the evolution of another set of random initial velocity gradient 
tensors. Every velocity gradient tensor in this set of data was normalized by the initial 
mean discriminant of the whole flow field of N total number of data points such that: 

Aii — Kj/D 1,€ > I>=T,D„/N, D^ +Q K (2.53) 

Ay is the initial velocity gradient tensor obtained from equation 2.52. The final 
velocity gradient tensor was obtained in the same way as described in the earlier 
section. However, each velocity gradient tensor was made dimensional by its own 
initial characteristic time scale t 0 (refer to equation 2.12) to obtain the dimensional 
velocity gradient tensor, Ay. The results are shown in figures 2.7 to 2.10. 

Figures 2.7(a-d) show the contour plots of the joint pdf of Q vs R at four differ- 
ent times. At initial time t = 0.0, the contour plot shows that the data points are 
distributed fairly evenly on both sides of the Q-axis. However, due to the random 
way in which the velocity gradient tensors were constructed, tensors with complex 
principal eigenvalues (in which case the data points will lie above the dividing curve) 
out number those tensors with all three real principal eigenvalues (which lie below the 
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Figure 2.3: Time evolution Q vs R. Aij normalized by local discriminant and local 
i? 1 / 3 so that all data points must lie on the line R = 1.0. (a) t=0.0. (b) t=5.0. 
(c) t=10.0. (d) t=20.0. 
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Figure 2.4: Time evolution of Q s vs R s . Aij normalized by local discriminant and 
local i? 1 / 3 . (a) t=0.0. (b) t=5.0. (c) t=10.0. (d) t=20.0. 
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Figure 2.5: Time evolution of — Q s vs Q w . Aij normalized by local discriminant and 
local i? 1 / 3 . (a) t=0.0. (b) t=5.0. (c) t=10.0. (d) t=20.0. 
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Figure 2.6: Time evolution of Q w vs (R s - R). A tj normalized by local discriminant 
and local R 1 ? 3 . (a) t=0.0. (b) t=5.0. (c) t=10.0. (d) t=20.0. 
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boundaries below the dividing curve). Figure 2.7(b) shows a very interesting result. 
There is an obvious shift of the data points toward the lower right quadrant. The 
overall structure is remarkably similar to the one obtained from the direct numerical 
simulation of a viscous , homogeneous isotropic flow shown in figure 4.8(d). The impli- 
cation is that even though the Restricted Euler model is based upon the assumptions 
that the effects of both the viscosity and the off-diagonal pressure may be neglected, 
the solution may still be applicable in certain regions of a real flow. However, the 
presence of the viscosity in real flow tends to prevent the flow field quantities from 
becoming singular. In this solution, the bulk of the data shifts slowly towards the 
lower right quadrant of the Q-R plot as time increases. 

The contour plots of the joint pdf between the second and third invariants of the 
symmetric tensor of this data set are shown in figure 2.8. At initial time t = 0.0, 
all the data points were evenly distributed below and within the boundaries of the 
dividing curve. As the velocity gradient tensor evolved from t = 0.0 to t = 0.2, 
an interesting result was observed. Most of the contour lines were relatively straight, 
extending from the upper left of the dividing curve towards the lower right. This same 
feature was also observed in DNS of a time-developing mixing layer[31], temporally 
evolving plane wake[30] and the homogeneous isotropic flow shown in figure 4.9. This 
linear characteristic of the contour lines persists until later times. 

Figure 2.9 shows the relationship between — Q s and Q w at different times. For most 
velocity gradient tensors, — Q s is comparable in magnitude to Q w initially. However, 
as the tensor evolves with the Restricted Euler model, — Q s begins to dominate Q w . 
This can be explained easily using the definition Q = Q s + Q w ■ Qs is negative definite 
while Q w is positive definite by definitions. As Q approaches — oo as dictated by the 
Restricted Euler model, Q s -4 Q while Q w -4 0, hence — Q s Q w . 

The sequence of evolution for Q w and ( R s — R) at four different times is shown 
in figure 2.10. At initial time t = 0.0, the data points formed a “tear-drop” shape 
structure, centered at the origin of the plot. As the solution evolves in time, the bulk 
of this structure slants towards the positive side of the ( R s — R) axis. This result is 
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consistent with the evolution equation for enstrophy density shown by Cantwell [4]: 

= 2 (Rs - R). (2.54) 

The next chapter will allude to the results presented here by way of comparison 
with a computation of a flow which is inviscid but includes the pressure. This will be 
followed in later chapters by comparisons with viscous simulations. 
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Figure 2.8: Time evolution of logarithmic contour plots of joint pdf of Q s vs R s . Aij 
is not normalized, (a) t=0.0. (b) t=0.2. (c) t=0.5. (d) t=1.0. 
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Figure 2.9: Time evolution of logarithmic contour plots of joint pdf of — Q s vs Q w . 
Aij not normalized, (a) t=0.0. (b) t=0.2. (c) t=0.5. (d) t=1.0. 
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Figure 2.10: Time evolution of logarithmic contour plots of joint pdf of Q w vs (R s —R). 
Aij is not normalized, (a) t=0.0. (b) t=0.2. (c) t=0.5. (d) t=1.0. 



Chapter 3 


Inviscid interaction of vortex tubes 


3.1 Introduction 

After observing the behavior of a random distribution evolving with the Restricted 
Euler model, it is of interest to determine if any region within a turbulent flow exhibits 
similar behavior. An inviscid calculation of interactions between two perturbed anti- 
parallel vortex tubes was therefore selected for this study. 

There are many debates on the existence of singularities in three-dimensional, 
incompressible Euler (inviscid) flow. Earlier inviscid simulations on two prescribed 
anti-parallel vortices by Pumir and Siggia [25] suggested that there is no singular- 
ity and the growth in peak vorticity at late times is at most exponential. Inviscid 
calculations with random initial conditions in a periodic box performed by Herring 
and Kerr [13] also appeared to grow at most at an exponential rate. Some viscous 
calculations [18, 20, 22] also suggested exponential growth when the trends of the 
calculations were extrapolated towards the limit of zero viscosity. Despite all these 
claims against the existence of singularities, Kerr [17] performed a simulation which 
provided evidence that supports the existence of a singularity. The simulation was of 
the interaction of perturbed anti-parallel vortex tubes using smooth initial profiles in 
a bounded domain. It was concluded by Kerr that the peak vorticity, the peak axial 
strain and the enstrophy density production scaled as ( t c — t) _1 while the enstrophy 


43 



44 


CHAPTER 3. INVISCID INTERACTION OF VORTEX TUBES 


IB- 


dividing plane 


symmetry plane 
Figure 3.1: The anti-parallel vortex tubes at initial condition. 


density grows logarithmically. t c was the critical time in the flow at which it was 
hypothesized to approach singularity. 

The objective of this thesis is not to address the debate on whether a singularity 
exists in Euler calculations, but to study the characteristics of this very interesting 
inviscid flow simulated by Kerr. The domain surrounding the peak vorticity in this 
flow was analyzed by classifying its local flow topologies. A coherent structure within 
this domain has been identified to exhibit similar characteristics predicted by the 
Restricted Euler model. 


3.2 Brief description of the flow 

Interaction of perturbed anti-parallel vortex tubes was simulated using DNS. Smooth 
initial profiles in a bounded domain with bounded initial vorticity were used to en- 
sure symmetry in the flow. A sketch of this flow at the initial condition is shown 
in figure 3.1. The symmetries imposed were represented by free-slip boundary con- 
ditions at the “dividing” plane between the vortices and the “symmetry” plane of 
the maximum perturbation of the vortices. Both planes are identified in the sketch. 
The image vortex tube, shown in dashed lines in the figure, was included to impose 
symmetry condition. 
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3.2.1 Numerical method 

Spectral methods that employed symmetries were employed in the simulation. The 
numerical method used Chebyshev polynomials which satisfy the free-slip boundary 
conditions at the dividing plane in the ^-direction, Fourier transforms in the periodic 
x-direction, as well as sine and cosine transforms in the y direction that satisfy free-slip 
boundary conditions at the symmetry plane. The mesh size used in the calculation 
was 512:256:192 in a domain of Air x 47r x 2it, which was concluded to have the best 
resolution properties. Global quantities such as enstrophy density and enstrophy 
density production were determined to have time dependencies consistent with the 
trends for peak vorticity and rate of strain. 

3.2.2 Initial condition 

Initial vorticity profiles used in the calculations with a high-wave number filter are 
similar to those used in Melander & Hussain [21] and Kerr & Hussain [18]. However, 
the use of the Chebyshev method in the vertical z direction imposed several limita- 
tions. Kerr overcame the difficulties by using a rather sophisticated way of defining 
a compact vortex core such that the unfiltered vortex core goes to zero smoothly at 
a given radius r, but is still Gaussian in the center. A high order exponential filter 
was then imposed to smooth out the rough edges. This filtering also expanded the 
vortex such that its edge was just above the dividing plane. 


3.3 Classification of local flow topologies 

In the results presented below, the domain of the flow analyzed is in a 27T x 27t x 7t box. 
This domain is only a small part of the whole flow, surrounding the peak vorticity 
in the flow. The initial peak vorticity w p \ 0 was set to 1.0 at the beginning of the 
simulation. This vorticity, with dimension of characteristic frequency, was used to 
non-dimensionalize time t obtained from the simulation: 

t = tx w p \ 0 . 


(3.1) 
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Data were analyzed from r = 12.5 when this domain of flow began to develop, 
until the latest time of r = 17.5 when this domain was concluded to exhibit self- 
similar characteristics. The non-dimensionalized critical time at which the domain 
was hypothesized to approach singularity is r c = 18.9 ±0.1. Significant resolution 
effects began to appear in the results after r = 17.0. However, the result at r = 
17.5 was included in the analysis to infer the scaling of the quantities that develop 
asymptotic behavior. In the discussion of results in the following sections, emphasis 
will be placed on the data set at time r = 17.0 when the convergence of all the strain 
components is strongest and the resolution most reliable. 

Figure 3.2 shows a contour plot of the iso-enstrophy density of the perturbation in 
the domain analyzed at time r = 17.0. This plot is presented from a perspective view 
looking along the stream-wise y-direction. The vertical z-axis has been stretched to 
five times its original value to show the development in the vertical direction more 
clearly. The horizontal x-axis is the dividing plane which separates the corresponding 
perturbation from the image vortex (not shown) directly below the dividing plane. 
Kerr labeled the region extending above the peak vorticity just behind the leading 
edge of the vortex as “head” while the vortex sheet extending behind the peak vor- 
ticity as “tail” . Both the “head” and “tail” were concluded to be vortex sheets, and 
remained so as the flow evolved to later times. The region of peak vorticity concen- 
trates into a region where the vortex sheet of the head and tail meet at a sharp angle 
at the lower left corner. 


Figure 3.3 shows the time evolution of the invariants Q and R from r = 12.5 
when the flow began to develop, until r = 17.5. The invariants “explode” from 
the origin at early time since the velocity gradient scales ~ l/(r c — r). Figure 3.4 
also exhibit similar observations where the invariants grow very rapidly from the 
origin at early time. This huge difference in the scale of invariants at different times 
makes comparison of invariants relatively difficult. A normalizing factor based on 
the instantaneous mean strain-rate was therefore chosen to normalize the velocity 
gradient: 


Aij = AC/Sl; Q = 


( (sjjSjjY 


1/2 


SijSij 


(3.2) 
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A\ 3 is the dimensional velocity gradient obtained from the simulation. 

Figure 3.5 shows the time evolution of the joint pdf between the second and third 
invariants of the normalized velocity gradient tensor. The contour plot of the joint 
pdf of the invariant approaches a self-similar shape. All the data points were observed 
to be at relatively low R values, almost hugging the Q-axis. No strong tendency of 
data points towards the right discriminant curve is observed. 

Two distinct structures are observed from the Q„-R s plots shown in figure 3.6. One 
of the structures consists of data points that hug the Q^-axis. The other structure is 
more complicated as it is actually a combination of different substructures. One of the 
substructures consists of data points that share the same ratio of a : /? : 7. Figure 3.7 
re-plots the same invariant plots with curves of fixed a : 0 : 7 ratio superimposed 
on the structures. Normalizing the eigenvalues by |/?|, the ratio of 7 : a ranges from 
-1.287 at r = 12.5 to -1.732 at r = 17.5. 

The other substructure identified consists of data points with the linear relation- 
ship: 

A 3 + XQ S + R S = 0. (3.3) 

A is any one of the three principal eigenvalues of the strain-rate tensor, S t] . This 
particular line osculates the discriminant curve at a single point and has the property 
that all the data points that lie on it have a common A. The identification of data 
points that lie on this line in both the invariant and physical spaces will be discussed 
in section 3.3.1. 

Figure 3.8 plots the relationship between — Q s and Q w of the velocity gradient 
tensor. Most of the data points in the plots have comparable local enstrophy den- 
sity and local dissipation, a direct result of data obtained from regions of high local 
enstrophy density due to the vortices and also high local strain due to the interac- 
tion/stretching between the two vortex tubes. Numerous data points in this flow have 
very high local dissipation but very low local enstrophy density. Another interesting 
observation from this figure is that data points with high local enstrophy density but 
very low local dissipation (eg. solid body rotation) are not identified in this flow. 

Figure 3.9 reveals very interesting plots showing the relationship between the 
local enstrophy density and the vortex stretching terms in the flow. The invariant 



3.3. CLASSIFICATION OF LOCAL FLOW TOPOLOGIES 


49 




Figure 3.3: Time evolution of logarithmic contour plots of joint pdf of Q vs R for 
un-normalized Aij. (a) r = 12.5. (b) r = 15.0. (c) r = 17.0. (d) r = 17.5. 
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Figure 3.4: Time evolution of logarithmic contour plots of joint pdf of Q s vs R s for 
un-normalized Aij. (a) r = 12.5. (b) r = 15.0. (c) r = 17.0. (d) r = 17.5. 
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plots show that all the data points are bounded within a fan-like structure, with no 
isolated data points outside the two distinct boundaries on both sides. Cantwell [4] 
concluded from the Restricted Euler model that there is a direct relationship between 
the evolution of local enstrophy density and the vortex stretching term: 

% = 2 (R, - R). (3.4) 

This relationship indicates that for data points evolving with the Restricted Euler 
model, there is an upper limit as to how fast the local enstrophy density of these 
data points can evolve in the flow. This limit is imposed by the vortex stretching 
term (R s — R). In addition, data points which satisfy the Restricted Euler model 
also approached an asymptotic state with relationship shown in equation 2.49 and 
figure 2.1. Therefore, data points that lie in the distinct right boundary of figure 3.9(c) 
were investigated further. These are the identified data points for further studies. 

3.3.1 Classification of local flow topologies for identified data 
points 

Figure 3.10 shows the various joint pdf for data points that are identified to lie on the 
right boundary of figure 3.9(c). The slope of Q w to ( R s — R ) for these data points is 
determined to be = 2.30. 

The Q-R and Q s -R s invariant plots are very interesting because all the data points 
identified can be best-fit onto a straight line on both plots. The straight line is of the 
form A 3 + A Q + it! = 0 in the Q-R plot and A 3 4- A Q s + R s = 0 in the Q s -R s plot. 
These lines osculate the respective discriminant curves at a single point. A is any one 
of the three principal eigenvalues of the velocity gradient tensor (in the case of Q-R 
plot) or the rate-of- strain tensor (in the case of Q s -R s plot). Therefore, all the data 
points that fall exactly on this line share the same eigenvalue. But which of the three 
principal eigenvalues ( a , (3 or 7) does A correspond to ? To answer this question, the 
three principal eigenvalues of both the Aij tensor (two of which may be complex) and 
the Sij tensor (all three of which must be real) of these identified data points were 
determined and sorted in descending order such that a > (3 > 7. These eigenvalues 
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(a) (b) 



(c) 


(d) 


Figure 3.6: Time evolution of logarithmic contour plots of joint pdf of Q s vs R s . 
(a) r = 12.5. (b) r = 15.0. (c) r = 17.0. (d) r = 17.5. 
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Figure 3.8: Time evolution of logarithmic contour plots of joint pdf of — Q s vs Q w . 
(a) r = 12.5. (b) r = 15.0. (c) r = 17.0. (d) r = 17.5. 



56 


CHAPTER 3. INVISCID INTERACTION OF VORTEX TUBES 



Figure 3.9: Time evolution of logarithmic contour plots of joint pdf of Q w vs (R s — R). 
(a) r = 12.5. (b) t = 15.0. (c) r = 17.0. (d) r = 17.5. 
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are plotted in figure 3.11. Looking first at the eigenvalues for the Sij tensor, there is a 
clear indication in figure 3.11(a) that the positive intermediate principal eigenvalue, 
0, is the eigenvalue that is shared by most of the data points. The mean value of 
0 was determined to be 0.40. Figure 3.11(b) plots the principal eigenvalues of the 
Aij tensor. At least one of the three eigenvalues is real, and that is the one shared 
by most of the data points. Data points with complex eigenvalues can usually be 
identified by those that have the same 0 and 7 in the plot. It is clear from this plot 
that the eigenvalue of the A^ tensor shared among the data points is the positive a. 
The mean value of a was also determined, and interestingly, it was also found to be 
0.40. Therefore, A corresponds to a in the Q-R plot and the positive 0 in the Q s -R s 
plot. The straight lines with A = 0.40 are superimposed in both figures 3.10(a) and 
(b) to show how well the data points fit onto the line. 

In the discussion of the behavior of the velocity gradient tensor in an incompress- 
ible flow where the acceleration gradient tensor Hij ^ 0, Cantwell[5] hypothesized the 
existence of an intermediate asymptotic state for the velocity gradient tensor A l j(t). 
Cantwell theorizes that after a fluid particle is set into motion by any flow, the particle 
settles into a state similar to the Restricted Euler model where its velocity gradient 
tensor satisfies the asymptotic form: 

MV » Kij[R(t)f 3 . (3.5) 

The angular momentum of this fluid particle changes relatively slowly under the action 
arising from the flow when Hjj / 0. Cantwell proceeded to propose the intermediate 
asymptotic state for Aij(t): 

Ai(i) = MyeS'W; ^ = -4«/(f). (3.6) 

For particles which satisfy this intermediate asymptotic model, a particular relation- 
ship that links the discriminant of H^, ,D|#.. to that of A^, D |^.. was derived: 

0|», i =DU, J (fi + /Q + / s ) 2 (3.7) 

The main consequence of this relationship is that the discriminants of and H v 
have the same sign, since the factor in brackets is squared. The eigenvalues of both 
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gradient tensors were therefore concluded to have the same character. The result 
that data points from the incompressible, Euler calculations were observed to have 
the same relationship as the factor in brackets lends support to the validity of this 
model. 

The — Qs-Qw and Q W -(R S — R) invariant plots of these data points also reveal very 
interesting results. These data points also exhibit characteristics of the asymptotic 
state of the Restricted Euler model; described in equations 2.47 and 2.49. Most of 
these data points appear to lie on a straight line, intercepting the — Q s axis at a 
finite positive value in the —Q s -Qw plot and the origin in the Q W -(R S — R) plot. The 
straight lines that best fit these data points are found to be: 

-Q s = 0.89Q W + 0.41; (3.8) 

Q w = 2.30(R S - R). (3.9) 

Similar data points were also identified in the same way at various times from 
r = 12.5 to t — 17.5. All the data points can be best-fit to a straight line with 
equation A 3 -I- A Q s + R s = 0. Figure 3.13 depicts these data points together with the 
straight line when A = 0.40. 

These same data points were also identified in the physical space as shown in 
figure 3.14, which is a contour plot of the local enstrophy density of the flow above 
the dividing plane. This figure illustrates the local enstrophy density level in the 
symmetry plane. The contour levels range from 0.445 to 4.45 at interval of 0.445 
for all four plots. Data points that are identified in the right boundary have been 
superimposed on these plots, marked by the symbol ‘+\ These data points formed a 
very coherent structure from the earlier times when the flow began to develop until 
the later times; always located near the peak vorticity and evolved with it. 

3.3.2 Classification of local flow topologies for data points 
with different Q W /(R S — R ) slopes 

In addition to data points that are identified on the right boundary of the Q w vs 
(R s — R) plot of figure 3.9(c) (with slope Q W /(R S — R) = 2.30), other data points 
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Data points 

(a) 



Data points 

(b) 


Figure 3 . 11 : Principal eigenvalues of identified data points with Q W /{R S — R) == 2.30. 
x:a, o:,3 and + 17 . (a) rate-of-strain tensor Sij. 0 — 0.40. (b) velocity gradient tensor 
Aij. a = 0.40. 
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(c) (d) 


Figure 3.14: Time evolution of the contour plots of local enstrophy density in the 
symmetry plane. Contour levels range from 0.445 to 4.45 at interval of 0.445. +: 
data points that lie on the right boundary of the Q w vs ( R s — R ) plots, (a) r = 12.5. 
(b) r = 15.0. (c) r = 17.0. (d) r = 17.5. 
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that lie on different Q W /(R S — R) slopes at r = 17.0 are also identified and studied. 
These data points are within ±15% of slopes ranging from 5.00, —8.00 to the left 
boundary in figure 3.9(c), with slope of —3.80. Due to the extremely small number of 
data points that lie on the left boundary, the slope of —4.40 was selected to replace 
—3.80 instead. Choosing a slope of —4.40 ensures that a reasonable number of data 
points are provided, and still include data points that lie on the left boundary. Data 
points that lie very close to the Q w axis are also of interest. These data points have 
very large Q w f (R s — R) values. 

Figures 3.15 and 3.16 depict the Q-R and Q s -R s plots of the selected data points. 
Again, the points tend to lie on a straight line of: A 3 + A Q s + R s = 0. The value of 
A was found to vary depending on the slope. 

An interesting question arises as to which principal eigenvalue does A correspond to 
when the slope selected is negative ? To answer that, the principal eigenvalues of both 
the Aij and SV, tensors for selected data points with slope of —4.40 were determined 
and plotted in the same way as in figure 3.11 shown earlier. Figure 3.17(a) plots out 
the eigenvalues of the SV? tensor. These data points do not show as strong a tendency 
to share an eigenvalue as those identified with slope of 2.30 in the earlier section. 
However, most of these data points have roughly equal value of the intermediate 
principal eigenvalue f3 , which is negative. The mean value of f3 was found to be 
—0.26. When the principal eigenvalues of the A^ tensors for these data were also 
determined, only a portion of the data points appear to have roughly equal value of 
7- 

Similar studies were done for data points that lie on other slopes. These data 
points share roughly equal value of the intermediate principal eigenvalue of the rate- 
of strain tensor. This common principal eigenvalue is positive when the slope is 
positive, negative otherwise. 

On the other hand, these data points do not have strong indication of having 
the same principal eigenvalue of the velocity gradient tensor when the slope selected 
deviates from the distinct right boundary with slope of 2.30. Among those data 
points that have roughly equal eigenvalue, the eigenvalue that these data points have 
in common is a when the slope is positive, 7 otherwise. 
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Table 3.1 gives the corresponding A value to the different slope selected. A is 
determined from the mean (3 of the data points on these slopes. The same A is also 
used in the equation of the straight line superimposed onto Q-R plots for these data 
points. 


Table 3.1: Values of A corresponding to selected slopes. 


slope 

2.30 

5.00 

OO 

-8.00 

-4.40 

A 

0.40 

0.20 

« 0.0 

-0.16 

-0.26 


Figure 3.18 shows that most of these points tend to have comparable local dissipa- 
tion and local enstrophy while figure 3.19 identified the data points that correspond 
to the various slopes. These data points were also identified in the physical space 
shown in figure 3.20. Data points with positive slopes tend to lie closer to the “tail” 
of the vortex contour while those of negative slopes tend to lie close to the “head” of 
the vortex contour. 

Results from subsections 3.3.1 and sec: others indicated that data points that lie 
on the right boundary of the Q w vs (R s - R) exhibit characteristics predicted by the 
Restricted Euler model. These data points have roughly the same intermediate prin- 
cipal eigenvalue of the rate of strain tensor, as well as the largest principal eigenvalue 
of the velocity gradient tensor. The relationship of — Q s and Q w for these data points 
is also of the same form as the asymptotic solution of the Restricted Euler model. In 
the physical space, these data points form a coherent structure in the flow, always 
very close to the peak vorticity, These results indicate that even though Restricted 
Euler model assumes the acceleration gradient tensor Hij = 0, its solution may still 
be applicable and useful in certain regions of a “real” flow. Similar results were not 
observed for data points that lie on other Q W /(R S — R) slopes. 
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(a) 



(b) 


Figure 3.17: Principal eigenvalues of data points with Q W /{R S - R) = —4.4. x:a, o :/? 
and +:7- (a) rate-of-strain tensor S tr (b) velocity gradient tensor Atj. 
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Figure 3.19: Q w vs (R s — R) for data points with various Q W /(R S —R) slopes, (a) 5.00. 
(b) oo. (c) -8.00. (d) -4.40. 
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(a) 


(b) 



(c) 


(d) 


Figure 3.20: Contour plots of local enstrophy density in the symmetry plane. Contour 
levels range from 0.445 to 4.45 at interval of 0.445. +: data points with various 
Qw/{R S ~ R) slopes, (a) 5.00. (b) oo. (c) -8.00. (d) -4.40. 
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3.3.3 Classification of local flow topologies for other interest- 
ing structures 

In addition to studying the local flow topologies for data points with various slopes of 
Q w and (R s — R) as well as locating them at the corresponding physical space, data 
points that lie within other interesting structures were also identified. 


Data points with low local enstrophy density 

Figure 3.21 depicts the invariant plots of data points with local enstrophy density 
less than or equal to five percent of the maximum enstrophy density in the flow at 
time t = 17.0. Referring to figure 3.8(c), these are the data points with relatively 
high local dissipation. Figures 3.21(a) and (b) are almost identical. This is because 
Q — Q s + Q w by definition, which implies that Q ~ Q s since these data points have 
very low Q w . These data points were identified in the low contours of the vorticity 
contour plot in figure 3.25(a). 


Data points with high local enstrophy density 

The invariants of data points with local enstrophy density of more than 70% of the 
maximum enstrophy in the flow were also examined. Both Q-R and Q s -R s plots gave 
early indication that all these data points might have the same 7/a ratio. However, 
figure 3.22(b) reveals that these data points do not correspond to any of the ratios 
superimposed onto the plot. The reason is because the tendency for data points to 
have the same 7/a ratio is usually observed only in regions of high local dissipation. 
In this flow, data points with high local enstrophy density do not always have high 
local dissipation, as observed in figure 3.22(c). Therefore, these data points do not fit 
onto any of the 7/a ratio curves. Figure 3.25(b) identified these data points in the 
physical space, which all lie within the vortex core with the highest contour level at 
the lower left corner where the “head” meets the “tail”. 
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Data points with high local dissipation 

Two distinct structures were formed by data points with local dissipation greater than 
70% of the maximum dissipation in the flow shown in the Q-R plot in figure 3.23(a). 
The first structure with mostly positive Q and negative R is very similar to the one 
shown in figure 3.23(a). In fact, these are the data points that also have high local 
enstrophy, a subset of those data points discussed in the preceding subsection. The 
other structure is formed by data points with high local dissipation but with low local 
enstrophy density. Two separate structures are also formed by these data points in 
figure 3.23(b). One of the structures hugs the Q s axis, while the other fits onto a 
curve of common 7 /a ratio shown in figure 3.7(c). These data points are indicated 
in the physical space shown in figure 3.25(c). 

Data points with comparable local enstrophy and dissipation 

Finally, data points of local enstrophy density within 15% of the local dissipation are 
identified in figure 3.24. These data points have very low Q, as shown in the Q-R 
plot. This is because Q s is always negative while Q w is always positive. Therefore, 
<5 for data points with exactly the same magnitude of Q s and Q w would be zero 
by definition. Figures 3.24(c) and (d) show that these points formed a very coherent 
structure in both plots while figure 3.25(d) identified these data points in the physical 


space. 







3.3. CLASSIFICATION OF LOCAL FLOW TOPOLOGIES 


75 



(a) 


5. Or 


-QstAjpU 


0.0 ! 1 i i i i i j l i i 

0.0 Q w (Ay) 5.0 



(b) 


5.0r 


Q w (Ajj ) L 


n n l i i i i i i r i i i 

-5.0 Rs(Aij)-R(Aij) 5.0 


(c) 


(d) 


Figure 3.22: Data points with Q w > 70% of Q w | maz . (a) Q vs R. (b) Q s vs R s . 
(c) -Q s vs Q w . (d) Q w vs ( R s - R). 
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(C) (d) 


Figure 3.25: Contour plots of enstrophy density in the symmetry plane. Contour levels 
range from 0.445 to 4.45 at interval of 0.445. +: data points with (a) Q w < 5%Q w \ max - 
(b) Q w > 70 %Q w \max. (c) IQs I > 70%|Qs|mox- (d) -Qs ~Qw 



Chapter 4 

Homogeneous isotropic flow 


A random field is homogeneous if its statistical properties are invariant under a trans- 
lation of the coordinate system. The one-point statistics of a homogeneous field are 
independent of position, and multiple-point statistical properties depend only on the 
separation of the sampling points. A random field is isotropic if all its statistical 
properties are invariant under a rotation of the coordinate system. The one-point 
statistics of such a field are independent of the orientation of the coordinate system. 
Turbulent fluid motion is usually homogeneous if it is isotropic. Although ideal ho- 
mogeneous and isotropic turbulence rarely exists, regions within turbulent flows often 
exhibit characteristics of homogeneous and isotropic flow. 

This chapter focuses on a direct numerical simulation of a homogeneous and 
isotropic flow in a periodic box. This relatively low Reynolds number flow simulation 
was performed on a Silicon Graphics Indy workstation. 


4.1 Approach 

The three dimensional Navier-Stokes equations for incompressible flow were solved nu- 
merically using a pseudo-spectral method. The flow was assumed to be homogeneous 
and isotropic, with periodic boundary conditions imposed on all sides of a cube. The 
length of each side of this cube is 2tt. The computational procedure involves trans- 
forming the velocity field from the physical space to Fourier space. Calculations of 
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solutions to the Navier-Stokes equation are performed in this spectral wave-vector 
space. The convolution sum that appears in the Fourier transformed Navier-Stokes 
equation is most efficiently calculated using a discrete form of the convolution the- 
orem and the Fast Fourier algorithm. The non-linear products are first calculated 
in physical space and then transformed back to the spectral wave-vector space. A 
straight forward application of the Fast Fourier transform introduces “aliasing” er- 
rors. This is because the Fast Fourier transform is a discrete Fourier transform with 
no distinction made between wave-vector components of modulo 2mk c ; where k c is 
the maximum wave- vector (Nyquist “frequency”) considered in the calculations and 
m is any integer. Energy of high wave-vector components unresolved in the simula- 
tion is thus “aliased” into components of lower wave numbers. One way of removing 
the aliasing errors more efficiently is to truncate the interactions outside a specified 
boundary in Fourier space. This method does not significantly affect the resolution 
because only the high-wavenumber tail of the spectrum, containing a small fraction 
of the total energy, is truncated. 


4.1.1 Derivation of governing equation 


Given a velocity field in physical space, u = u(x, t), the direct Fourier transform of 
u from physical space to spectral wave-space is: 


u(k, t) 


(2tt) 


J [exp[-ik • x]u(x, f)dx, 


(4.1) 


where u(k, t) denotes the transformed velocity field in spectral wave-space. To convert 
u(k, t) to physical space, apply the inverse Fourier transforms to u(k, t): 


u 


(x, t) = J [exp(ik • x]u(k, t)dk. 


(4.2) 


Both integrals involved in the direct and inverse transforms are carried out over an 
infinite three-dimensional domain. In subsequent sections, the symbol will be 
used to denote the transformed quantities in the Fourier wave-space. 

Consider the incompressible Navier-Stokes where the continuity and the momen- 
tum equations are, respectively: 


V • u = 0, 


(4.3) 
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— + u • Vu = -VP + i/V 2 u. (4.4) 

P = p/p is the kinematic pressure and v = pj p is the kinematic viscosity. 

The Fourier transform of the continuity equation is: 

V • u = 0 = k • u(k, t ) (4.5) 

Equation 4.5 reveals an important result: the Fourier transform of the velocity, 
u, is constrained to lie in a plane in wave-space, tt, perpendicular to k since the 
dot-product between k and u is zero at all times. 

Taking the Fourier transform of the momentum equation: 

- + vk 2 u(k, f) = — zkP(k, t) — u • Vu (4.6) 

or 

d 

(— + vk 2 ) u(k, t) = -ikP(k, t) - u • Vu (4.7) 

The terms on the left hand side of the above equation involve only the velocity, 
u, which evolves only in the wave-space plane tt. Therefore, all the terms on the 
right hand side must also evolve in the same wave-space plane. However, kP(k, t) is 
parallel to k and perpendicular to the plane -k. Hence, this quantity must be equal 
and opposite to the out-of-plane component of u • Vu. This constraint was used to 
remove the pressure term from the momentum equation. Since ikP(k,t) + u • Vu 
must lie in it - plane, this quantity can be rewritten as: 

k • (ikP 4- u • Vu) = 0, (4.8) 


or 


ik 2 P + k • u Vu = 0. 


Multiplying equation 4.9 through by i , an expression for P is obtained: 

ik -VL ■ Vu 


P = 


k 2 


(4.9) 


(4.10) 


Finally, the momentum equation can be rearranged to obtain the governing equa- 
tion for u: 
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The quantity u • Vu is not written as a convolution sum in wave-space because it 
is most efficient to carry out the product in physical space during the numerical im- 
plementation. This method of calculating the non-linear term is the pseudo-spectral 
method (also known as the Fast Fourier transform method). 

4.1.2 Numerical method 

The governing equation for Ui in equation 4.11 was simplified further so that a nu- 
merical algorithm could be implemented to solve it more efficiently. 

Treatment of the non-linear term 

The non-linear term u • Vu is the convolution sum, responsible for most of the algo- 
rithmic complexities of the numerical scheme. Using the expression TT to denote the 
direct Fourier transform and TT~ l for the inverse Fourier transform, the convolution 
sum could be replaced by fj m : 

(u^Vu) = TT{n-Vu)=ik m TT[TT-\u m )?T-\u j )) (4.12) 

fjrn 

Equation 4.11 becomes: 

U*(k, t) = £2 k — U • Vu = — ik m flm + (4-13) 

or 

dUl fr' — = ikmflm + ~ vk 2 u t {\s., t). (4.14) 

Treatment of the viscous term 

The numerical implementation of equation 4.14 could be simplified by removing the 
viscous term on the right hand side. This viscous term can be absorbed into the 
main source term using a very elegant technique employed by Rogallo[28]. Given an 
Ordinary Differential Equation (ODE) of the form: 


dill yO . 1 2 ~ 

— = G + vk u u 
at 


(4.15) 



4.1. APPROACH 


83 


this equation can be simplified by introducing an integrating factor, <f> 
that: 

d . , „ . dui o 

— (0u;) = <j>— + 4>uk 2 iii - <t>G. 

This same technique is used to simplify equation 4.14: 

f ^ u, = nk,y^f jm ~ u, m ] = Gi, 


dt 


{e ukH ui) = e 


vk 2 t 


i , Mkj , 
1'kmflm "b l „ f. 


k 2 


jm 


e uk2t such 


(4.16) 

(4.17) 

(4.18) 


Time marching scheme 

A second-order Adams-Bashforth time marching scheme was selected to discretize 
equation 4.18 forward in time. The Adams-Bashforth scheme converts an ODE from: 


dm „ 

~r = F - 
dt 


( 4 . 19 ) 


to 


At 


< +1 - uf = — [3 F n - F n ~ r . 


Applying the same scheme to equation 4.18 gives 


At, 


(^) n+1 ~ (#i) n = -y [3 {4>Gi) n - («^G/) n-1 ]. 


(4.20) 


(4.21) 


Substituting Q = e"* 2 ‘ into the above equation, 

e vk*{t n +Vt) £n+l _ e ukH » & n = ^[3 e ^t nG n _ e »k\t n -Vt) Q n - ( 4 . 22 ) 

2 

Dividing this equation by e" fc2tn • e" fc2vt , and rearranging: 


uf +1 = e-^^uf + ^[3e~ uk2At Gf - e-^^Gf- 1 }. 


(4.23) 


Gi — — Equation 4.23 is the evolution equation for velocity in wave 

space ui. 
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Evolution of the zero-wave number component 

An apparent singularity is observed in equation 4.23 at the zero wave number (where 
k = \Jkl + kj + kl = 0). This singularity can be removed by considering equation 4.6. 
At zero wave number, this equation becomes: 

= -ikP(0, t) - u^Vu(0, t) 

at 

= -ikjTT^T^iuj) ■ FT - 1 (in)] 

= -ikjfij( 0 ,t) 

= 0 

Therefore, the zero wave number component of the velocity does not evolve. This 
apparent singularity due to the zero wave number component can be avoided by 
considering only the non-zero wave number components in the evolution equation. 


4.1.3 Initial condition 


The velocity in the simulation was initialized to satisfy both the continuity and 
isotropy conditions. This initial condition was first used by Rogallo[28]. Given an 
initial energy spectrum, the continuity condition is satisfied by assuming velocity u: 

u = uiei = a(k)e[ + 0(k)e' 2 . (4.24) 


d is the computation vector basis and e- is any vector basis having e' 3 parallel to k. 
The complex components a and (3 in general are random in amplitude and phase, 
subject only to the constraint that: 


J <u • u* > dA(k) = E(k ) = J <aa* + (30* > dA(k). 


(4.25) 


E(k ) is the desired energy spectrum and the bracket <> denotes the inner product 
between the vectors. Rogallo selected a and 0 to be: 

a= (l§) ' e ' cos ^ ^(f!) e " si " 0 - (4 - 26) 


0i, 62 and <() are any random numbers uniformly distributed within the interval (0, 27 t). 
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To complete the formulation, the basis ' needs to be related to the computational 
basis e*j. This is achieved by using any basis subjected to the constraint: 

ke 3 ' = kyey + k 2 e 2 + k 3 e 3 = k. (4.27) 

Rotation about e 3 is absorbed into the random phase <f>. A particular basis is arbi- 
trarily chosen such that: 

ei ' • e 3 = 0. (4.28) 

This choice of basis leads to a solution for e) ': 


{k\ + ^) 1/2 e*! ' = k 2 e 3 - k x e 2 (4.29) 

k{ki + k 2 ) 1 c 2 = k\k 3 e\ + k 2 k 3 e 2 — (A; 2 + k 2 )e 3 . (4.30) 

Therefore, an initial velocity with an energy spectrum E(k ) that satisfies both the 
continuity and isotropy conditions is: 


( ®kk 2 + 0k\k 3 _ pk 2 k 3 — akki _ i P(k^ + k 2 ) x ^ 2 _ ^ 

V k(k\ + kl) 1/2 61 + k[k 2 + k%y7 2 e 2 + l e 3 J 


(4.31) 


In the simulation, the initial energy spectrum was chosen to be of the form “Top- 
hat”, as shown in figure 4.1. The energy for wave numbers k = ^Jk\ + k% + k% in 
the range between kj 3 and 2k c /Z were initialized to a uniform magnitude. Energy 
for wave numbers outside this range was set to 0.0. k c is the Nyquist wave number, 
which is half the number of grid points used in the simulation on each side of the 
cube. 


4.1.4 Accuracy and Stability analysis 

The numerical scheme used in this simulation is spectrally accurate in all the spatial 
directions, and second order accurate in time. 

The stability of the numerical scheme can be analyzed by applying the second 
order Adams-Bashforth time marching scheme to a model equation. This model 
equation is a one dimensional linear convection-diffusion equation: 


du du d 2 u 

dt a dx V dx 2 


(4.32) 
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m 



k 


Figure 4.1: Initial energy spectrum. 


Both a and v are constants in this equation. When 2tt = NVx, velocity u at a single 
Fourier component k can be expressed in a general form as: 

u = e at e 7 * ikx / NVx = e at ■ e ikx . (4.33) 


Substituting the expression of u into equation 4.32: 

^ = — (i /k 2 + iak)u = Lu. (4-34) 

at 

A second order Adams-Bashforth time marching scheme is used to discretize the 
above ODE into: 

u n+1 = u n + ^(3 Lu n - Lu n ~ l ). (4.35) 

Let a denotes e aVt . Since e at = (e“ vt ) n = a n , it is clear that the criterion for 
numerical stability is |cr| < 1. 

Substituting u n = a n e ikx into equation 4.35: 

<j 2 - (1 + \LVt)o + ^LVt = 0. (4.36) 

z z 


The roots for this equation are: 




1 + LVt + 7 L 2 Vt 2 
4 


(4.37) 
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Re (L Vt) 

Figure 4.2: Stability region for second order Adams-Bashforth marching scheme 
Using power series to expand the terms within the square-root: 

Jl + LVt + h^VP = 1 + l(LVt + h 2 W)-l(LVt + h 2 Vt 2 ) 2 + 

V 4 Z A 8 4 

^(LVi + ?I 2 Vt 2 ) 3 + ---. (4.38) 

The two roots of equation 4.36 are: 

oi = 1 + LVt + \l?Vt 2 + 0(Vt 3 ); 

£ 

0 2 = l -LVt - ^L 2 Vt 2 + 0(Vt 3 ). (4.39) 

The first root shows that Adams-Bashforth scheme is second order accurate in time. 
The second root 0 2 is spurious. The presence of spurious root is a characteristic of 
all multi-step time marching schemes. In the case of second order Adams-Bashforth, 
this spurious root is of less concern since 0 2 — > 0 as Vt — ^ 0. Figure 4.2 shows the 
stability region when the numerical scheme is stable. The stability contour crosses the 
Re(LVt) at -1 and tangent to the imaginary axis. Strictly, the method is unstable for 
purely imaginary L (when u = 0). However, this instability is very mild. In viscous 
flow, the scheme is stable. 
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4.1.5 Aliasing errors 

One of the complications arises in using the pseudo-spectral method to evaluate the 
convolution sum is the “aliasing” of high wave-number components into the lower 
ones. This is due to the fact that fast Fourier transform is incapable of distinguishing 
between wave-vector components of modulo 2mk c , where k c is the Nyquist wavenum- 
ber and m is any integer. In Fourier transform, each Fourier coefficient denotes the 
contribution of a specific wave number to the velocity represented. Aliasing refers to 
the error in some Fourier coefficients due to contamination by other Fourier coeffi- 
cients with wave numbers outside the range — k c < k < k c . Figure 4.3 illustrates this 
error when the number of grid points used in the simulation ( N = 8) cannot resolve 
the wave-vector components represented by the higher wave numbers outside the 
range of — k c < k < k c , (k c = Nf 2 = 4). Figure 4.3 shows that there are enough grid 
points to represent a sine wave of sin (30) but not sin (90). The information captured 
by the grid points made the sine wave sin(90) appears to be sin(0). 

In this simulation, the aliasing error is eliminated by the use of two sets of shifted 
grids as well as truncating Fourier wave components higher than a specific range. The 
details of this de-aliasing technique are described in Appendix B. 


4.2 Classification of local flow topologies of A{j ten- 
sor 


The characteristic time scale t 0 (eddy turn-over time) is used to normalize the dimen- 
sional time t obtained from the simulation: 


r = t/t 0 \ t 0 = q 2 0 /e 0 , 


(4.40) 


where Ql = \uiUi and e 0 = 2usijSij are the mean kinetic energy and the mean dissipa- 
tion at initial condition. Figure 4.4 plots the time history of the mean kinetic energy, 
the mean strain-rate and the mean Reynolds number based on Taylor microscale A. 
A is defined to be: 


A = 



Sij Sjj 


(4,41) 
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(a) 



(b) 


Figure 4.3: Aliasing error due to insufficient number of grid points used in Fourier 
transform. — : actual sine wave to be represented, ‘o’: Fourier coefficients represented 
by the grid points. “aliased” sine wave represented by the grid points, (a) sin(30). 
(b) sin(90). 
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Self-similar solutions obtained from dimensional analysis[3j in homogeneous shear 
flow indicates that length 6 and velocity u scale with time t as: 

<5 ~ t k ; u 0 ~ t k -\ (4.42) 

Therefore, 

Aij ~ u 0 /6 ~ t 1 . (4.43) 

Using this scaling relation for the length and velocity scales, the similarity solution 
dictates that the mean kinetic energy, mean strain-rate and the mean Reynolds num- 
ber should scale with time t in the following ways: 

q 2 ~ UiUi ~ t 2k ~ 2 , 

Aij Aij t , 

Re x ~u 0 - 5 ~ t 2k ~\ (4.44) 

These predicted scalings were superimposed on the plots with k — 0.31. 

The energy spectrum and other related spectra for the various velocity derivatives 
of this homogeneous isotropic flow are presented in Appendix C. 

Figure 4.5 shows a general view of a homogeneous isotropic flow in a periodic box 
obtained from the simulation. Figure 4.5(a) is a contour plot of iso-enstrophy density 
surfaces, with |cu| = 2.8. The characteristic length scales of the structures shown 
in this plot are relatively small, due to the relatively low Reynolds number in the 
simulation. This plot is compared to a contour plot of discriminant D in the same flow, 
shown in figure 4.5(b). Recall from the definitions of invariants that D = (27/4)R 2 + 
Q 3 , where Q and R are the second and third invariants of A i} respectively. Therefore, 
the same velocity gradient tensor Aij that was used to determine vorticity in the 
flow was also used to determine D in corresponding grid points. Figure 1.3 depicts 
various topologies in an incompressible flow, depending on D. For positive D , the 
local topology is either a stable- vortex/stretching or an unstable-vortex/ compressing. 
Figure 4.5(b) plots the iso-contour surfaces with D = 6.0. This plot reproduces 
structures obtained from contour plots of iso-enstrophy density surfaces. 

Figure 4.6 shows the time evolution of the invariant plots from r = 0.0 to r = 5.30. 
At initial time r = 0.0, the contour plot of this joint pdf is evenly distributed on both 
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(a) 


(b) 




(c) (d) 




(e) 


(0 


Figure 4.4: Time history of (a) mean kinetic energy, (b) log-log plot of mean ki- 
netic energy, (c) mean strain rate, (d) log-log plot of mean strain rate, (e) Taylor 
microscale Reynolds number, (f) log-log plot of Taylor microscale Reynolds number. 
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(a) 



Figure 4.5: Iso-contour surfaces for homogeneous isotropic flow at r = 5.30. (a) local 
enstrophy density |u;| = 2.8. (b) local discriminant D = 6.0. 
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sides of the Q-axis. As the flow evolves, the velocity gradient decays at a rate ss 1/r. 
Therefore, the invariants at later times are very small compared to the initial values. 
Figure 4.7 shows similar observations for the Q s -R s invariant plots. This rapid decay 
of invariants makes comparison of velocity gradient tensors at different times relatively 
difficult. To overcome this problem, the instantaneous mean strain-rate was used to 
normalize the velocity gradient tensor: 

= tJ , i /2’ (4-45) 

[SijSij ) 

A'ij is the dimensional velocity gradient tensor obtained from the simulation. 

Figure 4.8 shows the time evolution of the normalized Q-R invariant plots from r = 
0.0 to r = 5.30. As the flow’ evolves, the joint pdf of the invariants of normalized veloc- 
ity gradient approaches a self-similar shape. Data points with higher gradients tend to 
have local flow topologies of stable-focus/stretching and unstable-node/saddle/saddle. 
Similar trends have been observed in other turbulent flows [7, 30, 31]. Figure 4.9 shows 
the time evolution of the contour plots of the joint pdf of Q s versus R s . Data points of 
high gradients are again observed to lie close to the lower right quadrant of the Q s -R s 
plot. Since dissipation is directly related to Q s , these data points are concluded to 
reside in high local dissipation regions. Figure 4.10 indicates that local dissipation 
and local enstrophy density of this flow are comparable in magnitude. Figure 4.11 
depicts the tear-drop shape structure formed by the contour plots of the joint pdf 
of Q w versus (R s - R). This structure centers at the origin and tilts towards the 
positive ( R s — R) axis as the flow evolves. 


4.2.1 Effects of conditioning with local dissipation 

When data points are conditioned at higher local dissipation levels, there is a tendency 
for data points to move close to the right discriminant curve. These data points have 
strong preference for local flow topology of unstable-node/saddle/saddle as showm in 
figures 4.12 and 4.13. 



94 


CHAPTER 4. HOMOGENEOUS ISOTROPIC FLOW 




Figure 4.6: Time evolution of logarithmic contour plots of joint pdf of Q vs R for A!^ 
at (a) r = 0.0. (b) r = 1.32. (c) r = 2.65. (d) r = 5.30. 
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(a) 


(b) 



(c) 


(d) 


Figure 4.7: Time evolution of logarithmic contour plots of joint pdf of Q s vs R s for 
A' tj at (a) r = 0.0. (b) r = 1.32. (c) r = 2.65. (d) r = 5.30. 
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(C) 


(d) 


Figure 4.10: Time evolution of logarithmic contour plots of joint pdf of — Q s vs Q w 
for Aij at (a) r = 0.0. (b) r = 1.32. (c) r = 2.65. (d) r = 5.30. 



4.2. CLASSIFICATION OF LOCAL FLOW TOPOLOGIES OF A u TENSOR 99 


6.0r 


Q w (Ajj)k 


0.01 L 

- 6.0 


6.0r 


Qw^Ajj)k 


O.Ol L 

- 6.0 


6.0r 



6J0 


(a) 


(b) 


6.0r 



RgCAij^^Aij) 


i 


6.0 


Q w (Ajj)| 


O.Ol L 

- 6.0 



6.0 


(c) 


(d) 


Figure 4.11: Time evolution of logarithmic contour plots of joint pdf of Q w vs (R s -R) 
for A{j at (a) r = 0.0. (b) r = 1.32. (c) r = 2.65. (d) r = 5.30. 
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Figure 4.13: Logarithmic contour plots of joint pdf of Q s vs R s for A tj . r = 5.3. Data 
conditioned at various levels of maximum local dissipation, (a) 0% (All data points), 
(b) 25%. (c) 50%. (d) 75%. 
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4.2.2 Strain-rate distribution 

The three principal eigenvalues (strain-rates) of the symmetric tensor of Aij were 
determined and sorted in descending order such that: 

a > 0 > 7. (4.46) 

Each eigenvalue was normalized by the magnitude of the intermediate principal eigen- 
value, \0\ so that the intermediate principal eigenvalue is either 4-1 or -1. Figure 4.14 
shows the probability density functions of these normalized eigenvalues. The first plot 
includes all the data points while the second plot includes only data points with local 
dissipation higher than 25% of the maximum local dissipation in the flow. Most of 
the data points are observed to have positive 0. This observation is more pronounced 
with the conditioned data. The maxima for a and 7 are « 1.5 and —2.5 for the 
conditioned data. 

4.2.3 Vorticity-strain alignment 

Ashurst et a/.[l], Sondergaard[30] and Tsinober et al.[ 32] observed strong preferential 
alignment of the vorticity vector with the intermediate principal strain-rate direction. 
This alignment was stronger when data was conditioned at higher local dissipation 
regions. Figure 4.15 plots the probability density functions of the alignment angle 
between the vorticity vector and the three principal eigenvectors of the rate-of-strain 
tensor of this homogeneous isotropic flow. The vorticity vector cj, is related to the 
A^ tensor by: 

u \ = A32 — A23; 

U>2 = A13 — A31, 

U>3 — A 21 — Ai2- (4-47) 

The alignment angle, 0 is obtained from the dot product between this vorticity vector 
and each of the three principal eigenvectors of the rate-of-strain tensor corresponding 
to the principal eigenvalues a, 0 and 7. The vorticity vector is aligned exactly with 
the eigenvector when cos0 = 1.0(0 = 0) and non-aligned when cos0 = 0(0 = 90). 
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Figure 4.14: Probability density function of principal strain-rates normalized by |/?|. 
t = 5.30. —.a. (a) All data points, (b) |Q s (Ay)| > 25%\Q s {A ij )\ max . 
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Figure 4.15 shows that the vorticity vector tends to align with the principal eigen- 
vector corresponding to eigenvalue /?. This preferred alignment is more pronounced 
when data are conditioned with local dissipation greater than 25% of the maximum 
local dissipation in the flow. This same vorticity vector is observed to non-align with 
the principal eigenvector corresponding to 7, the most compressive strain-rate direc- 
tion. This non-alignment of the vorticity vector with the most compressive strain-rate 
direction was also observed by Sondergaard[30] in the plane wake simulation. 


4.3 Validation of the pressure Hessian tensor 


The pressure p at every grid point in the flow is determined using equation 4.10. 
Differentiating this pressure twice spatially gives the pressure Hessian tensor, P tJ : 


>, _ &P 

lJ dxidxj 


(4.48) 


Pij tensor is symmetrical by definition. Unlike the Aij tensor, the trace of this P tJ 
tensor is not zero, as indicated by equation 2.3 which gives the relationship between 
A^ tensor and the trace of the pressure Hessian tensor. This relationship is used to 
validate the numerical calculation of pressure by plotting Q(Aij) vs 0.5 P(Pij) (both 
equal to -%A ik A ki ) at r = 0 in figure 4.16. P(Pij) is the first invariant (trace) of Pij: 


\ nr.) A a2p 


2dx t dx,~ 2 AitAt ‘~ Q(Ail) ' 


(4.49) 


Notice that the data points do not all fall exactly on the Q{Aij) = 0.5 P{Pij) 
straight line. This is because Q(Aij) was obtained by multiplying the A^s in the 
physical space; whereas the pressure derivatives were determined in the Fourier space 
before being transformed to physical space. The wavenumbers included in the dis- 
crete Fourier transforms are limited within the range of ±N/ 2, where N is the number 
of grid points used in the simulation. Therefore, the calculated pressure derivatives 
contain only information from wavenumbers within ±iV/2. On the other hand, in- 
dividual A^ in physical space also contains wavenumbers within the range of ±N/2. 
However, the products of Ai k A k i are performed in the physical space, which increases 
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(a) 



(b) 


Figure 4.15: Probability density function of cosine of angle 6 between vorticity vector 
and principal strain-rate directions, r = 5.30. — : a. — : 0. • ■ •: 7. (a) All data 
points, (b) \Q t {Aij)\ > 25%\Q s (A,j)\ max . 
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the effective range of wavenumbers represented by the products. Therefore, A ik A k i 
contains information of higher wavenumbers outside the range ±N/2. Hence, there 
is slight difference between them. To compare 0.5 P{Pij) to Q(A ZJ ), the correlation 
coefficient between them was determined. This correlation coefficient was determined 
to be 0.965, with correlation coefficient of 1.0 being exactly the same. 


4.4 Classification of local flow topologies of Hij 

The characteristics of H^ tensor is of the most interest in this thesis. Equation 2.4 
shows that the behavior of velocity gradient tensor A io is governed by H v tensor, 
referred to as acceleration gradient tensor. 

To study the behavior of the A { j tensor, Cantwell [4] derived the analytical solution 
of Aij when H l3 = 0 and proposed an intermediate asymptotic model[5] for A tJ 
when H^ # 0. Ultimately, the questions that need to be answered are: What are 
the characteristics of ? How does it behave in a turbulent flow? Are there any 
similarities and differences between the local flow topologies of H i3 and A l} tensors? 
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There are three invariants for each H xj tensor: 

P(Hij) = H u = 0; 

Q(Hij) = ——H mn H nm ] 

H(Hij) ^HrnnHnkHkin. (4.50) 

This Hij tensor may be decomposed into a symmetric and an anti-symmetric 
tensor: 


Hij — Sij\h + Wij\ h ] 

SiiU = \(H (j + HjiY 

I h ~ — Hji). (4.51) 

The second invariants of the S lJ \ h and Wy | h tensors are represented by Q s {H i3 ) and 
Qw{Hij ) respectively. Unlike the invariants of A 13 tensor, the exact physical meanings 
of Q s (Hij) and Q w ( H l3 ) are still unclear. The principal eigenvalues of the symmetric 
SijU tensor will also be determined and referred to as the acceleration strain-rates (as 
opposed to the principal eigenvalues of the A i3 tensor which are known as the true 
strain-rates). Similarly, the equivalent vorticity vector will also be determined from 
Hij tensor using corresponding definition for the true vorticity vector obtained from 
the Ai 3 tensor: 

k'Vii = H32 — H23', 

Uh.2 — H\2 ~ H 2 l ] 

^hz = H$\ — H 13 . (4.52) 

Equation 2.5 gives the definition of H i3 from the pressure derivatives and the 
second derivative of Ai 3 . This same H tJ can also be obtained using equation 2.4. 
The results obtained from both equations were compared by plotting Q(Hij) versus 
R(Hij) in figure 4.17. The H l} tensor obtained from equation 2.5 is exact within 
the wave frequency range of ±N/ 2. All the calculations were performed using the 
Fourier coefficients before being transformed to physical space. On the other hand, 
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(a) (b) 


Figure 4.17: Logarithmic contour plots of joint pdf of Q vs R for H tJ obtained using 
(a) time derivatives and product of A ty (b) pressure cross-derivatives and second 
derivative of Aij . 


the Hij tensor obtained from equation 2.4 contained the the products of AikAkj ob- 
tained from the physical space. Therefore, these products contain information of 
higher wave numbers outside the range of ±iV/2, as explained in the earlier section. 
However, equation 2.4 requires the information of dAij/dt. This time derivative of 
A^ is determined using a second order accurate numerical method. Therefore, the 
accuracy of H^ tensor determined from equation 2.4 is limited by the accuracy of the 
numerical method used in getting the time derivatives. 

In the rest of this chapter, the Hij tensor is determined from equation 2.5. Calcu- 
lation of H^ using this equation is preferred because it requires only information at a 
single time. The H \j tensor is also normalized by the instantaneous mean strain-rate: 


Hij = 


a 

(SijSij)' 


(4.53) 


H[j is the dimensional gradient obtained from the simulation. 
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Figure 4.18 depicts the time sequence of Q(Hij) versus R(Hij) from r = 0.0 to 
r = 5.30. Almost all the data points are below the discriminant curve initially, 
indicating that tensor is symmetric with real principal eigenvalues. This H tJ 
tensor is completely dominated by the symmetrical pressure cross-derivatives when 
initial Reynolds number is large. As the flow evolves, the joint pdf of the invariants 
slowly transforms into the skewed tear-drop structure shown in the final plot. The 
contour plot of joint pdf of Q versus R does not approach any self-similar struc- 
ture. Instead, the invariants of this Hij tensor grow relatively larger as the flow 
decays. Data points with high gradients are observed to move towards the lower 
left quadrant in the Q-R plot. These data points have equivalent local flow topolo- 
gies of stable-node/saddle/saddle. These observations are opposite to the results 
obtained for Aij tensor shown in figure 4.8. The tendency for data points to lie 
on the lower left quadrant of the Q-R plot is also observed in the relatively higher 
Reynolds number simulation of an evolving planar wake (Chapter 5) and a Burgers 
vortex (Appendix D). 

Figure 4.19 shows the time evolution of the joint pdf of the second and third 
invariants of the acceleration rate-of-strain tensor. As the flow evolves, the joint 
pdf of the invariants skews toward the lower left quadrant. The magnitude of the 
gradients also grows relatively larger as the flow decays. The higher contour levels 
form straight lines extending towards the lower left quadrant of the Q s -R s plot. The 
equivalent local flow topologies for these data points are stable-node/saddle/saddle. 
The intermediate principal eigenvalues of this acceleration rate-of-strain tensor are 
negative for these data points. 

Figure 4.20 shows that Q s is much larger than Q w at earlier times when the H v 
tensor is symmetric. These two invariants become comparable in magnitude as the 
flow evolves. Figure 4.21 shows the contour plot of the joint pdf of Q w versus (R s -R). 
The contour plot “explodes” from the origin and forms the characteristic “tear-drop” 
shaped structure. This structure tilts towards the negative ( R s - R) axis at later 
times. 
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(a) (b) 



(c) 


(d) 


Figure 4.19: Time evolution of logarithmic contour plots of joint pdf of Q s vs R $ for 
Hij at (a) r = 0.0. (b) r = 1.32. (c) r = 2.65. (d) r = 5.30. 
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Figure 4.20: Time evolution of logarithmic contour plots of joint pdf of — Q s vs Q w 
for Hij at (a) r = 0.0. (b) r = 1.32. (c) r = 2.65. (d) r = 5.30. 
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Figure 4.21: Time evolution of logarithmic contour plots of joint pdf of Q , 
for Hij at (a) r = 0.0 (b) r = 1.32 (c) r = 2.65 (d) r = 5.30 
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Figure 4.22: Logarithmic contour plots of joint pdf of Q vs R for (a) —(d 2 p/ dxidxj — 
(S ij /3)d 2 p/dx k dx k ). (b) vd 2 A lj /dx k dx k . r = 5.30. 


4.4.1 Contribution from components of 


In addition to studying the time evolution of the invariants of Hij tensor, the rel- 
ative contribution from the components of H r] tensor was also analyzed. The H t] 
tensor consists of the trace-free pressure derivatives term and the third derivative of 
the Aij , expressed in equation 2.5. Figure 4.22 compares the Q-R invariant plots 
of the trace-free pressure cross-derivatives and the vd 2 Aij/dx k dx k tensors. The in- 
variant plot for the pressure cross-derivatives indicates that this tensor is symmetric 
(d 2 p/dxidxj = d 2 p/dxjdxi). Therefore, all data points lie below the discriminant 
curve. The invariant plots of Q s versus R s , — Q s versus Q w , and Q w versus (R s — R) 
for this pressure term are redundant since Q = Q s , R = R s and Q w = 0. On the 
other hand, the Q-R invariant plot for the vd 2 A tJ /dx k dx k tensor shows that this is 
the term that contributes all the “rotational” part of the H \j tensor. 
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4.4.2 Effects of conditioning with local dissipation 

How does the behavior of Hij tensor change when data points are conditioned with 
different dissipation levels ? Figures 4.23 and 4.24 show the contour plots of the 
joint pdf of Q versus R and Q s versus R s when data points are conditioned with 
local dissipation levels of 25%, 50% and 75% of the maximum local dissipation in 
the flow. Data points with higher local dissipation tend to move closer towards the 
left discriminant curve. These data points are observed to lie along this discriminant 
curve when conditioned at the highest local dissipation levels. 

4.4.3 Acceleration strain-rate distribution 

The three eigenvalues of the symmetric part of the H tj tensor were determined and 
sorted in descending order such that: 

a h >Ph> lh (4.54) 

These principal eigenvalues will be referred to as the acceleration strain-rates. The 
probability density functions for these acceleration strain-rates (each normalized by 
the magnitude of the intermediate acceleration strain-rate, 1 / 3 * 1 ) were determined 
in the way described for the strain-rate distribution of A i} . Figure 4.25 shows the 
probability density functions of the acceleration strain-rate distribution. When data 
points are conditioned with local dissipation greater than 25% of the maximum local 
dissipation in the flow, the intermediate acceleration principal eigenvalue is predomi- 
nantly negative. These data points have preferred local acceleration strain topology of 
stable-node/saddle/saddle. The distributions for a* and 7 * remain relatively broad, 
with small peaks at approximately 2.5 and —1.5 respectively. 

4.4.4 Equivalent vorticity - acceleration strain alignment 

An equivalent vorticity vector is defined from the components of the H i3 tensor using 
similar definitions for the true vorticity vector: 


Uhi — H32 — # 23 ; 
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Figure 4.23: Logarithmic contour plots of joint pdf of Q vs R for H i} . r = 5.3. Data 
conditioned at various levels of maximum local dissipation, (a) 0% (All data points), 
(b) 25%. (c) 50%. (d) 75%. 
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Figure 4.24: Logarithmic contour plots of joint pdf of Q s vs R s for H^. r = 5.3. 
Data conditioned at various levels of maximum local dissipation, (a) 0% (All data 
points), (b) 25%. (c) 50%. (d) 75%. 
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Figure 4.25: Probability density function of acceleration strain-rates normalized by 
\f} h \. r = 5.30. — : a h . P h . • • •: 7h- (a) All data points, (b) |Q s (Ay)| > 
25%| Qs(A{j) |mai- 
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U} h 2 = H 13 — H%i ; 

w/»3 = FI 2\ - H\2- (4.55) 

Figure 4.26 shows the cosine of the angle between the equivalent vorticity vector and 
the three acceleration strain-rate directions. Both plots indicate that there is prefer- 
ential alignment of the equivalent vorticity vector with the intermediate acceleration 
strain-rate. This result corresponds to the earlier observation found for the vorticity- 
strain alignment in the same flow. However, unlike the vorticity-strain alignment 
results, the equivalent vorticity tends to non-align with ot (Recall that the vorticity 
vector tends to non-align with the most compressive strain, 7 ). 
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(a) 



(b) 


Figure 4.26: Probability density function of cosine of angle 9 between equivalent 
vorticity vector and the acceleration strain-rate directions, r = 5.30. — : a^. — : 
j3 h . •••: 7 h . (a) All data points, (b) > 25%\Q a (Aij)\ ma x- 





Chapter 5 

Temporally evolving plane wake 


5.1 Introduction 


Sondergaard[30] performed a direct numerical simulation on a temporally evolving 
plane wake to study the effects of the initial conditions on the evolution of the plane 
wake; including the development of three-dimensionality, the mean flow, structure, 
growth rate and the mean properties of the far-wake turbulence. To overcome the 
limitations on the resolution due to restricted computer resources, Sondergaard em- 
ployed a temporal formulation, where the spatially evolving flow is approximated by 
a temporally evolving flow as it convects downstream. This formulation approximates 
the view of an observer convecting downstream with the flow structures of interest. 
Therefore, “time” in the temporal formulations becomes the measure of the level of 
development of the flow instead of “space” in the case of spatial formulation. By 
using temporal formulations, Sondergaard was able to achieve better resolution of 
the smaller scales. 

Certain assumptions are made in all temporal formulations. One of them is that 
the stream-wise rate of change is small compared to the scale of the structures being 
studied. Hence the mean flow may be approximated as being locally parallel. The 
other assumption made is that the computational domain can safely be limited to 
only a few large structures to achieve better resolution of the smaller scales. Doing 
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so confines the development of the large scale motions to a computational box which, 
at late times in the evolution, may only be a few wake widths long. 


5.2 Brief description of the flow 

The main intent of this work by Sondergaard was to examine the effect of the choice 
of the initial conditions on the development of the incompressible plane wake. All the 
flows were started from a laminar base flow. A small number of disturbance modes of 
specific wavelengths including fundamental, subharmonic and sub-subharmonic were 
then superimposed onto the base flow. These disturbance modes were chosen to be 
the most unstable eigenfunctions as predicted by linear theory. 


5.2.1 Numerical method 

Calculating free shear flows using spectral methods in turbulent flows which are peri- 
odic in two directions and have a non-periodic compact vorticity in the third direction 
posed many challenges. Different approaches have been used to approximate the in- 
finite direction and incorporating the spectral method at the same time, all with 
well documented draw backs. Sondergaard avoided the draw backs encountered by 
these previous methods with the use of an algorithm similar to the one presented 
by Corral & Jimenez[9] together with an unique way of dividing the non-periodic 
infinite direction into three domains. One of the domains contained all the vorticity 
in the flow while the other domains were vorticity free and extend from infinity to 
meet the vorticity-containing domain at the boundaries from above and below. The 
vorticity magnitude and vorticity gradient at the boundaries where the domains meet 
were forced to be zero. In this way, the vorticity-containing domain can be treated 
as fully periodic and can be solved with the standard pseudo-spectral technique. All 
the quantities in the simulation have been normalized by the initial wake half width 
and the initial free-stream velocity. A growing uniform grid was also used in the non- 
periodic infinite direction to maximize the resolution while keeping the vorticity at 
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the edges of the vorticity-containing domain small to satisfy the asymptotic matching 
condition. 

5.2.2 Initial condition 

The initial condition for this time developing plane wake consists of a Gaussian mean 
velocity profile in the stream-wise direction while the mean velocities in the other 
two directions were set to zero. In this way, the only non-zero mean vorticity profile 
is in the non-periodic infinite direction. Small two dimensional or three dimensional 
disturbances which are periodic in the stream-wise and span-wise directions were then 
added to the mean flow. These disturbances may be along any or all of the funda- 
mental, subharmonic or sub-subharmonic stream-wise wavelength; with disturbance 
phase angle <f> if desired. The disturbance phase angle is defined with respect to the 
two-dimensional fundamental disturbance, with 4> = 2n represents a physical shift of 
one fundamental wavelength. The amplitude of each disturbance eigenfunction was 
chosen such that the initial magnitude was found to be small enough for the initial 
wake growth to be within linear regime while at the same time large enough to allow 
the wake to become non-linear without excessive computational resources. 


5.3 Classification of local flow topologies of A tJ 

Flow simulations at two different Reynolds numbers were analyzed in this chapter. 
The first flow is a three dimensional simulation with a stream-wise disturbance added 
onto the subharmonic wavenumber in addition to the two dimensional disturbance 
along the fundamental mode. The Reynolds number for this flow based on the initial 
wake half width and the initial velocity defect is 346. This flow is identified as 
wk346 in the analysis. The other flow is another three dimensional simulation with a 
three dimensional disturbance along the subharmonic wavelength together with a two 
dimensional disturbance along the fundamental wavelength. The Reynolds number 
based on the same initial wake half width and initial defect velocity of this flow is 
2768. This flow' is identified as wk2768. Both flows have pairs of oblique waves with 
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equal and opposite span- wise wavenumbers oriented 60° to the span- wise direction. 
The choice of three-dimensional modes was motivated by the stability analysis which 
suggested that wave pairs at near 60° to the span-wise direction are the most unstable 
modes. 

The instantaneous mean strain-rate is used to normalize the velocity and acceler- 
ation gradients in the flow such that: 

U.( Sf (5.1) 

\ s ij s tj ) 

A[j and H[j are the un-normalized gradients obtained from the simulation. 

Figure 5.1 shows the iso-contour surfaces for enstrophy density and discriminant 
D = (27/4) R 2 + Q 3 for wk346 at time t = 102.7. The free-stream is flowing along 
the stream-wise x direction, as indicated in both plots. The vertical y-axis extends 
to infinity in both directions to capture all the vorticity in the domain calculated. 
The well defined span-wise rollers along the z-direction were depicted very nicely in 
both iso-contour surfaces. Figures 5.2 to 5.5 show the time evolution contour plots 
of the joint pdf of various invariants for wk346 at four different times. The basic 
characteristics of these figures are very similar to those obtained by Sondergaard for 
the same flow without the normalization shown in equation 5.1. The contour plot of 
the joint pdf “explodes” from the origin at early time and peaks at time t = 102.7 
before the magnitude of the gradients decreases as the flow approaches self-similar 
regime. The joint pdf maintains a familiar “tear-drop” structure with high gradient 
exhibiting a tendency towards the upper left quadrant of the Q-R plot (with local 
topology of stable focus-stretching) and the lower right quadrant (with local topology 
of unstable-node/saddle/saddle). 

The initial shape of the joint pdf of the Q s -R s invariant plot shows a very distinct 
structure determined by the initial condition. However, this structure grows towards 
the general shape observed in the later times as the flow evolves. Most data points 
with high dissipation have a tendency towards the lower right quadrant of the Q s -R s 
plot, almost hugging the right discriminant curve. 
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(b) 


Figure 5.1: Iso-contour surfaces for wk346 at t = 102.7. (a) local enstrophy density 
|o>| = 2.8. (b) local discriminant D = 5.0. 
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The joint pdf of the -Q s -Qw plot shows that data points have comparable local 
dissipation and local enstrophy density initially. However, the local enstrophy density 
began to dominate as the flow evolved and developed well defined span-wise rollers 
with relatively weak stream-wise stretching between the rollers. Once the rollers were 
established, the three dimensional disturbances grew rapidly and the stream-wise 
structures became more intense, under the influence of the straining field. There- 
fore, the magnitude of the local dissipation and the local enstrophy density becomes 
comparable again at the latest time shown. 

The joint pdf of the Q w -R s -R invariant plots show that data points with the high- 
est local enstrophy density usually occur at regions with positive vortex-stretching 
terms. These structures are very similar to those observed in the homogeneous 
isotropic simulation. 

Figures 5.6 to 5.9 show the time evolution of the joint pdf of various invariants 
for flow wk2768 at two different times; after the flow has fully developed and is 
approaching the asymptotic regime. The magnitude of the highest gradients and the 
fraction of the flow containing the high gradient motions are larger in this higher 
Reynolds number flow. Otherwise, the higher Reynolds number has very little effect 
on the general behavior and the overall shape of the joint pdf. The development cycle 
of exploding from the origin at initial condition and decaying as the flow approached 
self-similar regime was concluded by Sondergaard to be independent of Reynolds 
number. 


5.4 Classification of local flow topologies of 

The tensor for the wake flows is determined from equation 2.4, repeated here for 
reference purposes: 

% + + A,„A tj - A^AnM = Hi, (5.2) 

dtoxk 3 

Figures 5.10 to 5.13 show the joint pdf of the invariant plots of the Hij tensor for 
wk346 at two different times. The flow is fully developed and begins to decay. The Q- 
R plots show that most data points lie below the discriminant curve, indicating that 
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(a) (b) 



(c) 


(d) 


Figure 5.2: Time evolution of logarithmic contour plots of joint pdf of Q vs R for 
wk346, Aij. (a) t=22.8. (b) t=52.8. (c) t=102.7. (d) t=204.8. 
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Figure 5.3: Time evolution of logarithmic contour plots of joint pdf of Q s vs R s for 
wk346, Aij. (a) t=22.8. (b) t=52.8. (c) t=102.7. (d) t=204.8 
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(c) (d) 


Figure 5.4: Time evolution of logarithmic contour plots of the joint pdf of -Q s vs 
Q w for wk346, Aij. (a) t=22.8. (b) t=52.8. (c) t=102.7. (d) t=204.8. 
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Figure 5.5: Time evolution of logarithmic contour plots of joint pdf of Q w vs (R s - R) 
for wk346, Aij. (a) t=22.8. (b) t=52.8. (c) t=102.7. (d) t=204.8. 
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(a) 



(b) 


Figure 5.6: Time evolution of logarithmic contour plots of joint pdf of Q vs R for 
wk2768, Aij. (a) t=99.8. (b) t=194.6. 
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(a) 



(b) 


Figure 5.7: Time evolution of logarithmic contour plots of the joint pdf of Q s vs R s 
for wk2768, Aij. (a) t=99.8. (b) t=194.6. 
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(a) 



(b) 


Figure 5.8: Time evolution of logarithmic contour plots of joint pdf of -Q s vs Q w for 
wk2768, Aij. (a) t=99.8. (b) t=194.6. 
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Figure 5.9: Time evolution of logarithmic contour plots of joint pdf of Q w vs (R s - R) 
for wk2768, Aij. (a) t=99.8. (b) t=194.6. 
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these data points have symmetrical tensors. This is because when the Reynolds 
number is relatively high, the viscous effect becomes relatively less significant when 
compared to the effects of the pressure Hessian. Since the pressure Hessian tensor 
is symmetric, the Hij tensor tends to be symmetric as well. Therefore, most of the 
data points were observed to lie below the discriminant curve. At later time, the 
magnitude of the gradients are much smaller with most data points close to the 
origin. When examined closely (looking at the magnified plot), the data points in 
the higher contour levels are observed to move towards the lower left quadrant of the 
Q-R plot, exhibiting similar trend observed for the homogeneous isotropic flow. 

The joint pdf of the Q s -R s plots also exhibit trend of moving towards the lower left 
quadrant. The higher contour levels are almost linear and parallel, showing strong 
preference for topology of stable-node/saddle/saddle. Data points with high Q s (Hij) 
are not necessary the same data points that are found in the high local dissipation 
region, indicated by high Q t (Aij). The effect of local dissipation on the behavior of 
the tensor will be illustrated in the later subsections. Figure 5.12 shows that data 
points in wake flows have significantly higher equivalent local dissipation compare 
to the equivalent local enstrophy density. The structure formed by the joint pdf of 
Q W -(R S — R) is also similar to the “tear-drop” shape. 

5.4.1 Effects of Reynolds number 

Figures 5.14 to 5.17 show the invariant plots of the higher Reynolds number 
simulation wk2768. Data for this simulation were selected at close to the evolution 
times for wk346 for comparison purposes. The shapes of the joint pdf of various 
invariants remain very similar to those of wk346. The effect of higher Reynolds 
number increases the magnitude of the highest gradients, as reflected by the increase 
in axis scales. 

5.4.2 Effects of conditioning with local dissipation 

The behavior of the H tj tensor for data points that reside in different local dissipa- 
tion regions are shown in figures 5.18 to 5.25 for both flows. These data points are 
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Figure 5.10: Time evolution of logarithmic contour plots of joint pdf of Q vs R for 
wk346, Hij. (a) t=102.7. (b) t=204.8. (c) same as (b), with magnified axis scales. 
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(c) 


Figure 5.11: Time evolution of logarithmic contour plots of joint pdf of Q s vs R s 
wk346, H t j. (a) t=102.7. (b) t=204.8. (c) same as (b), with magnified axis scales. 
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Figure 5.12: Time evolution of logarithmic contour plots of joint pdf of — Q s vs Q w 
for wk346, Hij. (a) t=102.7. (b) t=204.8. 
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Figure 5.13: Time evolution of logarithmic contour plots of joint pdf of Q w vs (R s — R) 
for wk346, Hij. (a) t=102.7. (b) t=204.8. 
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(a) 



(b) 


Figure 5.14: Time evolution of logarithmic contour plots of the joint pdf of Q vs R 
for wk2768, Hij. (a) t=99.8. (b) t=194.6. 
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(a) 



(b) 


Figure 5.15: Time evolution of logarithmic contour plots of joint pdf of Q s vs R s for 
wk2768, Hij. (a) t=99.8. (b) t=194.6. 
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Figure 5.16: Time evolution of logarithmic contour plots of joint pdf of — Q s vs Q w 
for wk2768, Hij. (a) t=99.8. (b) t=194.6. 
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conditioned at different levels of local dissipation. The levels vary from 0% of the 
maximum local dissipation (all data points are included) to 75% of the maximum 
local dissipation in the flow. As the conditioning levels increases, the magnitude of 
the invariants of Hij are observed to decrease very rapidly. This observation indi- 
cates that the magnitude of the H i0 tensor gets smaller when the conditioning level 
of the local dissipation increases. This conclusion is arrived from the results in fig- 
ure 5.10(a), which shown that most Hij tensor is symmetric since most of the data 
points lie below the discriminant curve. From the definition of Q(Hij), when H l} 
tensor is symmetric, then 


Q(Hij) ~ H ik H ki ~ H ik H ik ~ \H i3 \ 2 . (5.3) 

When data is conditioned at higher local dissipation levels, Q(Hij) is very small, as 
observed in figures 5.18 to 5.25. Hence, magnitude of H i} gets smaller when data is 
conditioned at higher local dissipation levels. 

The result that the magnitude of H i3 tensor is small at high local dissipation 
regions is very interesting. This result indicates that in regions of high local dissipation 
(i.e. high velocity gradients), the magnitude of Hij can be small. In these regions, 
the analytical solution of the Restricted Euler model may be applicable. 


5.4.3 Effects of conditioning with local enstrophy density 

The behavior of the H l3 tensor for data points with different local enstrophy density 
are shown in figures 5.26 to 5.33. As the conditioning levels increases from 0% to 
75% of the maximum local enstrophy density, the magnitudes of the invariants of 
H^ are not observed to get smaller accordingly. Instead, data points with high local 
enstrophy density levels correspond to high gradients in all the plots shown. This 
result indicated that the magnitude of H X) is small only in the high local dissipation 
regions, but not necessary in the high local enstrophy density regions. 
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Figure 5.18: Logarithmic contour plots of joint pdf of Q vs R for wk346, Hij. t=102.7. 
Data conditioned at various levels of maximum local dissipation, (a) 0% (All data 
points), (b) 25%. (c) 50%. (d) 75%. 
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(c) (d) 


Figure 5.19: Logarithmic contour plots of joint pdf of Q s vs R s for wk346, H iy 
t=102.7. Data conditioned at various levels of maximum local dissipation, (a) 0% 
(All data points), (b) 25%. (c) 50%. (d) 75%. 
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(c) 


(d) 


Figure 5.20: Logarithmic contour plots of joint pdf of Q vs R for wk346, Hij. t=204.8. 
Data conditioned at various levels of maximum local dissipation, (a) 0% (All data 
points), (b) 25%. (c) 50%. (d) 75%. 
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(a) (b) 



(c) 


(d) 


Figure 5.21: Logarithmic contour plots of joint pdf of Q s vs R s for wk346, Hij. 
t=204.8. Data conditioned at various levels of maximum local dissipation, (a) 0% 
(All data points), (b) 25%. (c) 50%. (d) 75%. 
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Figure 5.22: Logarithmic contour plots of joint pdf of Q vs R for wk2768, H l3 . t=99.8. 
Data conditioned at various levels of maximum local dissipation, (a) 0% (All data 
points), (b) 25%. (c) 50%. (d) 75%. 
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(a) 


(b) 



(c) 


(d) 


Figure 5.23: Logarithmic contour plots of joint pdf of Q s vs R s for wk2768, t=99.8. 
Data conditioned at various levels of maximum local dissipation, (a) 0% (All data 
points), (b) 25%. (c) 50%. (d) 75%. 
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(c) 


(d) 


Figure 5.24: Logarithmic contour plots of joint pdf of Q vs R for wk2768, H tj . 
t=194.6. Data conditioned at various levels of maximum local dissipation, (a) 0% 
(All data points), (b) 25%. (c) 50%. (d) 75%. 



152 


CHAPTER 5. TEMPORALLY EVOLVING PLANE WAKE 



(a) 


(b) 



(c) 


(d) 


Figure 5.25: Logarithmic contour plots of joint pdf of Q s vs R s for wk2768, Hij. 
t=194.6. Data conditioned at various levels of maximum local dissipation, (a) 0% 
(All data points), (b) 25%. (c) 50%. (d) 75%. 
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(c) (d) 


Figure 5.26: Logarithmic contour plots of joint pdf of Q vs R for wk346, H iy t=102.7. 
Data conditioned at various levels of maximum local enstrophy density, (a) 0% (All 
data points), (b) 25%. (c) 50%. (d) 75%. 
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Figure 5.27: Logarithmic contour plots of joint pdf of Q s vs R s for wk346, H iy 
t=102.7. Data conditioned at various levels of maximum local enstrophy density, 
(a) 0% (All data points), (b) 25%. (c) 50%. (d) 75%. 
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Figure 5.28: Logarithmic contour plots of joint pdf of Q vs R for wk346, t=204.8. 

Data conditioned at various levels of maximum local enstrophy density, (a) 0% (All 
data points), (b) 25%. (c) 50%. (d) 75%. 
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(a) (b) 



(c) (d) 


Figure 5.29: Logarithmic contour plots of joint pdf of Q s vs R s for wk346, Hij. 
t=204.8. Data conditioned at various levels of maximum local enstrophy density, 
(a) 0% (All data points), (b) 25%. (c) 50%. (d) 75%. 
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Figure 5.30: Logarithmic contour plots of joint pdf of Q vs R for wk2768, Hij. t=99.8. 
Data conditioned at various levels of maximum local enstrophy density, (a) 0% (All 
data points), (b) 25%. (c) 50%. (d) 75%. 
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(a) 


(b) 



(c) 


(d) 


Figure 5.31: Logarithmic contour plots of joint pdf of Q s vs R s for wk2768, H iy 
t=99.8. Data conditioned at various levels of maximum local enstrophy density, 
(a) 0% (All data points), (b) 25%. (c) 50%. (d) 75%. 
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Figure 5.32: Logarithmic contour plots of joint pdf of Q vs R for wk2768, Hy. 
t=194.6. Data conditioned at various levels of maximum local local enstrophy density, 
(a) 0% (All data points), (b) 25%. (c) 50%. (d) 75%. 
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(a) 


(b) 



(c) 


(d) 


Figure 5.33: Logarithmic contour plots of joint pdf of Q s vs R s for wk2768, Hij. 
t=194.6. Data conditioned at various levels of maximum local (a) 0% (All data 
points), (b) 25%. (c) 50%. (d) 75%. 
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5.4.4 Acceleration strain-rate distribution 

It was observed by Sondergaard that when the strain-rates of the velocity gradient 
tensor were sorted in descending order a > (3 > 7; the ratio of these strain-rates was 
found to be a : (3 : 7 ~ 1.5 : 1.0 : —2.5 for data points conditioned in the high local 
dissipation regions. Positive /3 indicates that the local rate-of-strain topology is of 
the type unstable-node/saddle/saddle. This distributions is insensitive to Reynolds 
number and holds for all the three-dimensional wakes simulated by Sondergaard. 
Defining the principal eigenvalues of the symmetric part of the H \j tensor as the 
acceleration strain- rates, and sort them in descending order such that oth > Ph > Jh, 
the acceleration strain-rates distribution was also studied in a similar fashion. 

The probability density functions of these principal eigenvalues for wk346 and 
wk2768 are shown in figures 5.34 to 5.37. All the eigenvalues have been normalized 
by the magnitude of the intermediate principal eigenvalue \ j3 h \ so that the normalized 
intermediate principal eigenvalue can only take on values of ±1.0 

Two different plots are shown in each figure. The first plot includes all the data 
points in the flow while the other include points conditioned with local dissipation 
higher than 25% of the maximum local dissipation in the flow. More negative inter- 
mediate principal eigenvalue Ph than the positive ones are observed in the distribution 
functions, indicating that the local flow topology of stable-node/saddle/saddle is pre- 
ferred to the unstable-node/saddle/saddle. The distribution for ah and 7* are almost 
symmetrical on both sides of the distribution functions. No significant peaks are ob- 
served for the broad ah and 7 h distributions. Similar results are observed for flows at 
different times and Reynolds numbers. 

5.4.5 Equivalent vorticity-acceleration strain alignment 

Sondergaard reported preferential alignment of the vorticity vector with the inter- 
mediate rate-of-strain tensor direction and non-aligned with the most compressive 
rate-of-strain direction by studying the alignment probability density function of the 
alignment angle between the vorticity vector with the three principal strain-rates 
direction. 
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(a) 



(b) 


Figure 5.34: Probability density function of the principal acceleration strain-rates 
normalized by |/3ft| for wk346, t=102.7. — : a^. — : /3ft. • 7ft. (a) All data points, 
(b) | Q t (Aij)\ > 25%\Q s {A tj )\ max 
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Figure 5.35: Probability density function of the principal acceleration strain-rates 
normalized by \fa\ for wk346, t=204.8. — : ct*. — : fa. ■ ■ •: 7a, . (a) All data points. 
(b)\Q s (A ij )\>25%\Q s (A ij )\ max 
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(a) 



(b) 


Figure 5.37: Probability density function of the principal acceleration strain-rates 
normalized by |/?J for wk2768, t=194.6. — : a*. — : Qh. • • •: 7 h- (a) All data points, 
(b) |«.(4,)| > 25%|Q,(^)U I 
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To compare the results between the velocity gradient tensor and the acceleration 
gradient tensor, similar studies were done for the H{j tensor by studying the align- 
ment angle between the equivalent vorticity vector (as defined in equation 4.55) with 
the acceleration strain-rates >/?/,> 7 / 1 ). Figures 5.38 to 5.41 plot the probability 
density functions of the cosine of the alignment angle 6 between the equivalent vor- 
ticity with the three acceleration rate-of-strain directions of the Hij tensor for both 
wk346 and wk2768. 9 is any one of the three alignment angles between the equiv- 
alent vorticity and the acceleration rate-of-strain direction. The density functions 
do not show any significant difference between a h and 7 h in all the plots. However, 
there is a strong preference for the equivalent vorticity vector to align with the in- 
termediate rate-of-strain, fih direction. This tendency becomes more pronounced (as 
indicated by the increased density functions levels) when data are conditioned with 
local dissipation of more than 25% of the maximum local dissipation in the flows. 
This observation appears to be independent of Reynolds number. 
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(a) 



(b) 


Figure 5.38: Probability density function of cosine of angle 6 between the equivalent 
vorticity vector and the equivalent principal strain-rate directions for wk346, t=102.7. 
— : a h . 8 h . • • •: 7 h . (a) All data points, (b) |Q s (A i: ,)| > 25%\Q s (A ij )\ max 
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(a) 



(b) 


Figure 5.39: Probability density function of cosine of angle 0 between the equivalent 
vorticity vector and the equivalent principal strain-rate directions for wk346, t=204.8. 
— : a h . — : (5 h . ■ • y h . (a) All data points, (b) |Q,(j4y)l > 25%|Q s (Ay)| mai 
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(a) 



(b) 


Figure 5.40: Probability density function of cosine of angle 9 between the equivalent 
vorticity vector and the equivalent principal strain-rate directions for wk2768, t=99.8. 
— : a h - 3 h - ■■ 7 h- (a) All data points, (b) > 25%\Q s (A ij )\ max 
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(a) 



(b) 


Figure 5.41: Probability density function of cosine of angle 6 between the equiva- 
lent vorticity vector and the equivalent principal strain-rate directions for wk2768, 
t=194.6. — : a h . ---: Q h . •••: 7 h - (a) All data points, (b) \Q s (A i:j )\ > 

25%\Q s (Aij)\ max 





Chapter 6 
Conclusions 


The Lagrangian behavior of the velocity gradient tensor was studied using results 
from three direct numerical simulation of turbulent flows: an inviscid calculation of 
interactions between two vortex tubes, a homogeneous isotropic flow and a temporally 
evolving planar wake. Analysis of these flows revealed strong preference for the inter- 
mediate principal eigenvalue of the rate-of strain tensor to be positive, with preferred 
local strain topology of unstable-node/saddle/saddle. There is also a preference for 
the vorticity vector to align with this intermediate principal eigenvector. These prefer- 
ences were found to be stronger when data were conditioned at higher local dissipation 
regions. All these observations support earlier findings of other direct numerical sim- 
ulated flows, including homogeneous isotropic flow and homogeneous shear flow[l], 
mixing layer[31], channel flow[2] and turbulent boundary layer[32]. These flows are 
very different in terms of large-scale motions and Reynolds numbers, indicating that 
these observations may be universal characteristics that are flow-independent. 

A model of the evolution of A tJ based on the solution of a Restricted Euler equa- 
tion contained asymptotic behavior consistent with these observations. This solution 
was compared with results obtained from an incompressible Euler calculation of two 
interacting vortex tubes. A coherent structure within this inviscid calculation was 
observed to exhibit characteristics predicted by the Restricted Euler model. This 
coherent structure evolved with the flow as the calculation approached singularity, 
always near the peak vorticity in the flow. Invariants of data points within this 
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structure have linear relationships of the form: 

4- A a Q + R = 0, (6-1) 

A^ + A S Q S + R s = 0. (6.2) 

A a was concluded to be the most positive principal eigenvalue of the velocity gradient 
tensor while A s was the intermediate principal eigenvalue of the rate-of-strain tensor. 
Both A a and A s were found to be 0.40. The invariants of the velocity gradient tensor 
for data within this structure have relationships similar to those predicted by the 
Restricted Euler model. 

The tensor H^, which is the anisotropic part of the acceleration gradient tensor, 
was studied in two different viscous flows: a relatively low Reynolds number decaying 
homogeneous isotropic flow in a periodic box, and a temporally evolving plane wake 
at two different Reynolds numbers. The magnitude of the invariants of the Hij tensor 
of the homogeneous isotropic flow increased as the flow decayed. The intermediate 
acceleration strain-rate in this flow inclined to be negative, with a local strain topology 
of stable-node/saddle/saddle. Furthermore, the equivalent vorticity vector preferred 
to align with this intermediate acceleration strain-rate direction. 

Behavior of acceleration gradient tensor of a temporally evolving planar wake 
at two different Reynolds numbers was also studied. The intensity of the highest 
gradients increased when the Reynolds number of this flow was higher. Otherwise, the 
shapes of the contour plots of the joint probability density function of the invariants 
of H^ remained very similar. The H^ tensor was found to be nearly symmetric in this 
moderately high Reynolds number simulation. The intermediate acceleration strain- 
rate also inclined to be negative, with local strain topology stable- node/saddle/saddle. 
A strong preferential alignment between the equivalent vorticity and the intermediate 
acceleration strain-rate direction was also observed. The most interesting result from 
the study of this flow was that the magnitude of Hij was very small when data were 
conditioned at higher local dissipation regions. This result was not observed for the 
relatively low Reynolds number simulation of homogeneous isotropic flow. All the 
results obtained for the wake simulation were nearly the same at the two Reynolds 
numbers studied. 



Appendix A 


Velocity gradient tensor with 
random components 

A velocity gradient tensor Ay- needs to be constructed such that its components 
are obtained from a random number generator with mean E(x) = 0.0 and variance 
cf 2 = 1.0. To make this velocity gradient tensor more “realistic”, the tensor needs to 
be forced to satisfy the continuity constraint for incompressible flow. Furthermore, in 
the case of isotropic flow, the volume integral of the second invariant Q must be zero 
due to the balance of the pressure gradients acting on the surface of a control volume. 
Therefore, the constructed Ay would have to satisfy the condition that E(Q) = 0.0 
as well. 

To satisfy both conditions, the velocity gradient tensor Ay is constructed from 
the sum of a symmetric tensor 5y and an anti-symmetric tensor JTy such that Ay = 
Sy+Wy. The components of Sy and W XJ tensors are obtained from a random number 
generator, which generates random numbers x with mean E(x) = 0.0 and variance 
cr 2 = 1.0. Forcing the trace of the S XJ tensor to zero, the components of Sy are: 

S u = x l - + x 2 + x 3 ), 

S22 =x 2 - ^(xi + x 2 + x 3 ), 

S33 = x 3 — - (xi 4 - x 2 + x 3 ), 
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5 12 = S21 = X4, 

513 = 531 = £ 5 , 

*5*23 = S 32 = ^6- (A.l) 

Similarly, the components of the anti-symmetric tensor W i3 are: 


W n — W 22 — W33 = 0.0, 


W n = -W 21 = x 7 , 

W n = -W 31 = x s , 

W23 = —W 32 = Xg. (A. 2) 


xi,x 2 , ■ ■ . ,xg are nine different random numbers obtained from the random number 
generator. The Aij tensor constructed this way satisfies the continuity constraint. 

The second condition that the volume integral of Q must be zero can be satisfied 
by forcing the expected value of Q to zero, i.e. E(Q) = 0.0. Noting that 5^ tensor is 
symmetric and W tJ tensor is anti-symmetric, 

Q 2 A-imA m i — — S im Si m + ^WifnWirn (A. 3) 


Put 


Then 


Now, 


E(Q) = ~E(S im S im ) + \E(W im W im ) = 0.0. 
E(S im S im ) - E(W im W im ) = 0.0. 


(A.4) 

(A.5) 


E (SimSim) — E(S xl + 5x2 + *5x3 + S 21 + S 22 + S% 3 + S 31 + S 32 + S 33 ) 
= E(S 2 n ) + E(S 2 2 ) + E(S \ s ) + E(S 2 2l ) + . . . + E(S 2 3 ). 


For the case when E(x) = 0.0, 


o 2 = E(x 2 ) - E(x) 2 = E(x 2 ), 


(A.6) 
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and, for independent variables x, 

E{xiXj) = E(xi)E(xj) = 0.0; » ^ j. (A. 7) 

Therefore, from E(x f) = E(x\) = = <r 2 = 1.0, and E{x x x 2 ) = = 

E(x 2 x 3 ) = 0.0, then 

E{S im S im ) = 8a 2 . (A. 8) 

To make E(Q) = 0.0, the variance of every component of Wij tensor needs to be 
modified by a factor a, such that E{Q) = E(S im S im - W im W iTn ) = 0.0. Therefore, 

8cr 2 — 6o 2 a 2 = 0.0 => a 2 = (A. 9) 

o 

Finally, the velocity gradient tensor Aij with random number components which 
satisfies both the continuity and E(Q) = 0.0 constraints is obtained by adding S lJ 
and tensors together, such that: 

*$11 = x ,\ — -(£1 "b x 2 + £3)1 

S22 — %2 — 2 > ^ Xl *b X2 "b x 3 )i 
S33 = X 3 - + X 2 + X 3 ), 


S\2 — S 2 l = X4, 

S\3 = S31 = £5, 

£23 = 532 = X 6 - (A. 10 ) 


and, 


W u — W 22 — W 33 — 0.0, 


Wi 2 = — W 2 \ = \J 4 / 3 x 7 , 
IT 13 = — W 31 = \J 4/3x g , 
IT23 = — W32 = 4/3x9. 


(A. 11) 
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Appendix B 

De-aliasing technique 


Consider a three-dimensional convolution sum: 

M k ) = Y «(p)«(q); p + q = k. (B.l) 

|p|,|q|<fe c 

k c is the Nyquist “frequency” . This alias-free sum can be evaluated using only two 
sets of grids as well as truncating wave-vector modes outside a specified range, as 
described by Orszag[23]. Two sets of grids x(j) and x(j) are selected such that: 


x Ci) = r'0'i»j2.j3); o < < 2k c , m = 1 , 2 , 3 . 

K c 

x(j) = ^0'i + i,j 2 + ^,i3 + ^); 0 < i m < 2k c , m = 1,2,3. (B.2) 

The velocities represented on these two sets of grid points are: 

“(j) = Y «(k)exp[fk-x(j)], (B.3) 

|k|<fc c 

«(j)= Y £( k )exp[ik-x(j)], (B.4) 

\k\<k c 

v(j)= Y t)(k) exp[zk • x(j)], (B.5) 

|k|<fce 

Hi)= Y fi(k) exp[ik • x(j)]. (B.6) 

|k|<fcc 

The transforms of both iTl y and u • v are: 

™( k ) = Y «-5exp[-*fc-£(j)]; (B.7) 

0<j<2A: c 
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Figure B.l: Region of retained Fourier modes defined by D. 

tu( k ) = 7^3 u-vexp[-ik-x(j)}. (B.8) 

[ZK C ) 0 <j<2fc c 

Therefore, 

^[u)(k) + u)(k)] = rD(k) + ^u)(k + 2fc c n); (|k| < k c ). (B.9) 

The term involving n represents the spurious aliasing errors. Patterson [24] has shown 
that only two of the three components are independently ±1 while the remaining one 
is zero. Consequently, equation B.9 involves aliasing only to the extent of wave- vectors 
that are “doubly aliased” . 

To remove the doubly-aliased terms, Orszag[23] considered a region D such that 
u( k) = v(k) = 0 for k outside the region: 

D= {\k a \<k c a = 1,2, 3 ( B . 10 ) 

| \k Q ±k 0 \ < *§*• oi,f3= 1,2,3, a^/3 

This inequality defines the region D bounded by an 18-sided polyhedron, as shown 
in the perspective sketch of figure B.l. If k is outside the region D, then the doubly- 
aliased terms in equation B.9 must be identically zero since there are no nonzero 


179 


components of w(p), t)(q) that can contribute to tD(k+2fc c n) such that k+2/fc c n = p+q 
with k, p, q 6 D. Therefore, if the spectral truncation is made so that the only 
excitable Fourier modes lie within the region D, then three-dimensional alias-free 
convolution sums are computable exactly on just the two sets of grids x(j) and x (j). 
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Appendix C 

Energy and velocity derivatives 
spectra 


The energy spectrum for the homogeneous isotropic flow was determined and shown 
in figure C.l(a). This energy spectrum is plotted against the magnitude of the three 
dimensional wave number vector, k. Both the energy spectrum E(k) and the wave 
number k are defined as: 

E(k) = Ul (k ) 2 + u 2 (k ) 2 + u 3 (k) 2 , (C.l) 

and 

k = + &2 -b ^3 - (C. 2 ) 

ki, k 2 and k 3 are the wave numbers in the three spatial directions. The Kolmogorov 
k ~°/ 3 scaling was superimposed onto the same plot, for comparison. Only a very 
small portion of the energy spectrum exhibits the k~ h/z scaling, indicating that this 
low Reynolds number simulation is not close to the “turbulent” flow. However, the 
spectrum shows that the grid size used in the simulation is sufficient to resolve all the 
wave numbers. 

In calculating the velocity gradient tensor, Ay = dui/dxj, there is a concern that 
the grid size used in the simulation may not be able to resolve this velocity gradient 
at higher wave numbers. To validate the resolutions, the spectrum of Ay, denoted as 
Ea(^)i was determined in a similar way as the energy spectrum with E^ik) defined 
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as: 

Ea{^) — An {k ) 2 + A 22 (k ) 2 + A 33 (k) 2 . (C.3) 

The EA(k ) spectrum is plotted in figure C.l(b). This spectrum shows that the grid 
size used in the simulation is able to resolve the higher wavenumbers of the velocity 
derivatives. 

Higher derivatives of the velocity fields were required in calculating the Hij tensor 
from either equation 2.4 or 2.5. To get a “feel” for the resolution of these higher 
velocity derivatives, the spectra for dAijfdx m and d 2 Aij / dx m dx m were determined 
in a similar way as the energy spectrum. These spectra axe defined as: 




(C.4) 


EaA(k) = f-i -(k ) 2 + + o— 

dxmdxm Oxjndxjn ox m ox ri 


m = 1,2,3. (C.5) 


As the derivatives get higher, the higher wave numbers become more dominant since 
these derivative were obtained spectrally by multiplying the wave numbers to the 


velocity field. At higher wave numbers, the “energy” of these derivatives become 
much more significant compared to those obtained for the energy spectrum, as seen 
in the spectra shown in figure C.2. The abrupt cut-off of the spectra at the maximum 
k was due to the range of the wave numbers used in the simulation to remove aliasing 


errors, as explained in Appendix B. 
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(a) 



(b) 


Figure C.2: Spectra for the higher derivatives of the velocity gradient tensor for 
homogeneous isotropic flow at r = 5.30. t indicates the maximum wave number, k max , 
resolved in the simulation, (a) dAij/dxm spectrum, (b) d 2 A lJ /dx m dx m spectrum. 





Appendix D 

Acceleration gradient tensor of 
Burgers vortex 


Assuming a steady Burgers vortex with velocity field Uj(x, y, z): 


Ui - ~x - Re y( 1 - e r2 )/r 2 . 

(D.l) 

u 2 = ~^y + Re x(l - e~ r 2 )/r 2 . 

(D.2) 

u 3 = z. 

(D. 3 ) 


x, y and z are the three Cartesian coordinates with z being the vertical axis and 
r = yfx 1 + y 2 . Re = T /&itv is the Reynolds number where T and 1/ are the circulation 
and the viscosity of the flow respectively. 


D.l The velocity gradient tensor A tJ 


From the velocity field Uj(x,y,z), the velocity gradient tensor Ay = dui/dx-j was 
obtained by taking the corresponding derivatives of ui, U2 and U3 with x, y and 2: 


Aii — — + Re 


2 xy 


(l-e' r! -rV rJ ) 


(D. 4 ) 


A 12 



(1 — e -r2 ) + Re 




(D. 5 ) 
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A 2l = ^(l-e- r2 )-Re 


2x 2 


(1 -e~ r -r 2 e~ r )\ . 


. 1 o 

A 22 — Re 


2 xy 


L r 

A33 = 1 . 0 . 


(1 - e~ r2 - r 2 e _r2 )J 


^13 — ^23 — ^31 = ^32 = 0 . 0 . 


(D.6) 

(D.7) 

(D.8) 

(D.9) 


Notice that the Aij tensor is not symmetrical since A i2 / A 21 . Equations 1.16 to 
1.18 give the definitions of the various invariants of the tensor for incompressible 
flow. Since the first invariant P is zero for all incompressible flow, the second and 
third invariants were determined to be: 


Q(A,i) = -J - fle 2 [Jj(l - e -r * - rV' 2 ) 2 - e^]. (D.10) 

R(A„) = -i + JSe 2 [i(l - - r 2 e~' 2 ) 2 - e- 2 ' 2 ]. (D ll) 


D.2 The acceleration gradient tensor, Hr, 

The definition of the acceleration gradient tensor, Hij was derived in equation 2.4. In 
the case of a steady Burgers vortex, the expression for the tensor can be simplified 
to: 

Hi, = u k ^ + A it A kj - (A^A^f. (D.12) 

The various components of the Hy tensor were determined with the help of a sym- 
bolic mathematical software, MAPLE[6]. With the various derivatives involved, the 
expressions of the components of the tensor are very complex. However, employ- 
ing the fact that the Burgers vortex is an axisymmetric flow, these expressions could 
be simplified by studying only the case along any radial r-axis. Assuming x = r and 
y = 0, the various components of the H t] tensor are: 

H u = —j^(—28Re 2 +56Re 2 e~ r2 +Z2Re 2 r 2 e~ r2 -28Re 2 e~ 2r2 -32Re 2 r 2 e~ 2ri +3r 4 )/r 4 . 

(D.13) 
(D.14) 


H 12 = Re e" r \ 
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H 21 — 6 r (— 1 + 2r 2 ). (D.15) 

#22 = - ^ (20 Re 2 - 40 Re 2 e~ r 2 - 16i?e Ve^ 2 + 20#e 2 e" 2r2 + l§Re 2 r 2 e~ 2r2 + 3 r 4 ) /r 4 . 

(D.16) 

#33 = — 1 (4#e 2 — 8Re 2 e~ r7 — 8Re 2 r 2 e~ r7 +4Re 2 e~ 2r2 +8Re 2 r 2 e~ 2r2 — 3r 4 )/r 4 . (D.17) 

#13 = #23 =: #31 = #32 = 0.0. (D.18) 

Notice that the #y tensor is again not symmetrical since ff 12 ^ H 2 \. The invariants 
of have also been calculated using MAPLE. The first invariant P(#y) is again 
zero. The expressions for the second and third invariants are: 


Q{Hij) 


on 2 —2r 2 2 3 o2 _ 2r 2 Re 2 (l - e- r2 )e~ r2 

16 r 2 

1/2C 2 - 16/3#e 4 (l - e~ r2 ) 2 e~ 2r2 
r 4 

28 #e 4 (l -e- r2 ) 3 e- r2 _ 13 fle 4 (l - e~ r2 ) 4 
3 r 6 3 r 8 


(D.19) 


#(#^) = Re 2 e~ 2r7 r 2 - l -Re 2 € + ^Re 4 e~ 3r2 (l - e~ r2 ) + 
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20 


-2 4 
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(D.20) 


where 


C\ = Re 4 { 1 - e- r2 )V 2r \ (D.21) 

C 2 = Re 2 { l-e~ r2 ) 2 . (D.22) 
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There is an apparent singularity at r = 0 since both invariants Q and R involve 
terms with denominator containing r raises to various powers. This apparent singu- 
larity was removed using L’Hospital’s rule, giving: 

Q(«,)| r=0 = + iife 2 - ifie 4 . (D.23) 

Wii) l-o = -55 - - Afle 6 . (D.24) 

D.3 Behavior of Q{Hij) and R(Hij) with r 

Figure D.l(a) shows the variations of Q(Hij)/Re 4 with r for Reynolds numbers of 
2, 5, 10, 20 and 50. Q(Hij) is normalized by Re 4 so that a self-similar solution for 
Q(Hij)/Re 4 is obtained for large Reynolds numbers (Re > 50). A self-similar solution 
is also obtained for R(Hij)/Re G when the Reynolds number is large, as shown in 
figure D.l(b). 

Figure D.2 plots Q(Hij)/Re 4 against R(H t} ) / Re 6 , with the discriminant curve 
D = (27/4) it* 2 + Q 3 = 0.0 superimposed onto the plot. The self-similar solution for 
large Reynolds number begins on a point that lies exactly on the discriminant curve 
with Q(Hij) / Re 4 = —0.333 and R(Hij) / Re 6 = —0.074 at r = 0.0. As r increases, 
this solution osculates the other discriminant curve at Q(Hij)/Re 4 = —0.117 and 
R(Hij) j Re 6 = 0.015 before approaching the origin. 










Q(H ;j )/Re 
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Figure D.2: Q(Hij) / Re 4 versus R(Hij)/Re 6 . •••:Re = 2.- ■ :Re = 5. :Re = 10, 

— :Re = 20. — :Re = 50. 
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